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ABSTRACT 


A  Trimble  Advanced  Navigation  Sensor  (TANS)  (^uadrex  Global  Positioning  Sys¬ 
tem  (GPS)  receiver  processing  unit  and  three  antenna/preamplifier  assemblies  were 
flown  aboard  Space  Shuttle  Discovery,  Sl^-Sl,  as  part  of  DTO  700-6,  GPS  On-orbit 
Demonstration  (GOOD).  The  experiment  was  designed  to  quantify  advantages  and  iden¬ 
tify  potential  problem  areas  ftv  Space  Shuttle  GPS  operations  using  a  low  cost,  comn^r- 
cial,  space  configured,  GPS  receiver.  GPS  data,  including  position,  velocity,  time, 
health,  and  status  information  were  recorded  during  the  mission.  Following  the  mission, 
a  reference  trajectory  was  generated  by  NASA  Jdinson  Space  Center  through  post-pro- 
cessing  of  the  Orbiter’s  on  board  navigation  state.  The  recorded  GPS  data  has  been  ana¬ 
lyzed  and  compared  to  the  reference  trajectt^  to  evaluate  the  navigational  performance 
of  the  receiver.  Additionally,  post-flight  filtering  of  the  GPS  data  has  been  performed  in 
order  to  determine  whether  a  significant  increase  in  performance  may  be  obtained 
through  filtering. 
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I.  INTRODUCTION 


This  thesis  investigated  the  performance  of  a  low  cost,  commercial,  space  config¬ 
ured,  Global  Positioning  System  (GPS)  receiver^  flown  aboard  Space  Shuttle  Discovery, 
STS-51,  as  part  of  DTO  0700-6,^  GPS  On-orbit  Demonstration  (GOOD).  The  DTO  was 
sponsored  and  funded  by  NASA  Johnson  Space  Center,  and  developed  with  the  suf^xnt 
of  NASA  contractors,  and  students,  including  the  authors,  from  the  Naval  Postgraduate 
School.  Data  recorded  during  the  mission  was  analyzed  to  evaluate  navigational  perfor¬ 
mance  of  the  receiver.  Additionally,  post-flight  filtering  of  this  data  was  performed  in 
order  to  determine  whether  a  significant  increase  in  performance  could  be  obtained 
through  filtering. 

A.  BACKGROUND 

1.  Space  Shuttle  Orbiter  Baseline  Navigation  System 

A  wide  variety  of  equipment  is  employed  in  die  Orbiter’s  baseline  navigation 
system.  Ail  navigation  sensor  information  is  supplied  to  a  six-state  suboptimal  Kalman 
filter,  which  provides  the  navigation  functions  with  three  position  and  three  velocity 
states.  The  three  position  states  are  the  coordinates  specifying  the  Obiter’s  position  vec¬ 
tor  in  the  Aries-mean-of-19S0  (M50)  (Cartesian  coordinate  system.^  Likewise,  the  three 
velocity  states  define  the  Orbiter’s  velocity  vector  in  the  MSO  system. 


^Trimble  Advanced  Navigation  Sensor  (TANS)  (juadrex  GPS  Receiver  Processor  Unit. 
^Detailed  Test  Objective  Number  0700-6. 

^Tbe  MSO  system  is  defined  in  NASA  Technical  Memorandum  X-S81S3,  October  1974. 
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During  the  ascent  phase  of  a  mission,  the  Inertial  Measurement  Unit  (IMU) 
is  the  primary  sensor,  providing  attitude  and  acceleration  data  to  the  Kalman  Alter.  This 
data  is  augmented  by  ground  based  C-band  radar-tracking  information  uplinked  over  an 
S-band  communication  link.  During  the  on-orbit  coasting  phase  of  a  mission,  the  IMU 
provides  attitude  data,  and  acceleration  data  from  Orbital  Maneuvering  System  (OMS) 
burns.  Accelerations  falling  below  the  IMU  threshold  arise  from  opposing  Reaction 
Control  System  (RCS)  thrusters  that  do  not  form  a  perfect  couple  (vernier  effect),  and 
from  external  venting  of  gasses  and  waste  products.  These  unaccounted  for  accelerations 
result  in  steadily  increasing  navigtuional  error. 

While  on-orbit,  the  ground  continues  to  track  the  Orbiter  using  ground-based 
C-band  radar.  Two-way  Doppler  tone  ranging  over  S-band  and  Ku-band  communication 
links,  either  direct  or  via  the  Tracking  and  Data  Relay  Satellite  (TDRS)  system,  provides 
additional  tracking  capability.  When  die  Orbiter’s  navigational  state  is  observed  to  devi¬ 
ate  from  the  ground  based  tracking  trajectory  by  a  pre-defmed,  mission-dependent 
amount,  a  new  state  vector  consisting  of  three  position  states,  three  velocity  states,  and  a 
time  tag  is  uplinked  by  Mission  Control.  At  various  times  during  a  mission  (prior  to  ren¬ 
dezvous  and  deorbit  bum),  an  IMU  alignment  may  be  performed  using  an  on  board  star- 
tracker  to  correct  attitude  error  caused  by  gyro  drift. 

During  the  re-entry  phase  of  a  mission,  IMU  data  is  augmented  with  drag 
modeling  data  from  250,000  feet  down.  As  the  Orbiter  passes  through  the  ionosphere, 
all  radio-navigation  and  communication  signals  are  blacked  out  for  a  period  of  time. 
Upon  exit  from  blackout,  contact  is  first  established  by  C-band  tracking  radar.  A  state 
vector  can  be  uplinked  as  soon  as  S-band  communication  is  regained.  Subsequently,  the 
Orbiter  can  receive  L-band  Tactical  Air  Navigation  (TACAN)  station  signals  and  begin 
area  navigation,  combining  TACAN  range  and  bearing  data  with  barometric  (30,000- 
2,500  feet)  and  radar  (2,500  feet  down)  altimeter  measurements.  Final  approach  and 
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landing  are  accomplished  with  a  microwave  scan  beam  landing  system  (MSBLS),  begin¬ 
ning  normally  at  10,000  feet,  10  nautical  miles  downrange  from  touchdown. 

2.  Space  Shuttle  Orbiter  ^S  Navigation  System 

Pursuant  to  a  smdy  contract  commissioned  by  NASA,  Rockwell  Internation¬ 
al’s  Space  Systems  group  conducted  a  design  and  integration  smdy  of  a  GPS-based  pri¬ 
mary  navigation  system  for  the  Orbiter  in  the  late  1970’s.  (Van  Leeuwen  et  al,  1979,  pp. 
118-135)  The  smdy  demonstrated  that  the  use  of  on  board  satellite  GPS  receivers  for 
precise  orbit  determination  was  clearly  feasible,  and  expected  the  improved  navigation 
capabilities  to  yield  significant  (^rational  benefits.  The  smdy  concluded  the  GPS-based 
system  to  be  a  technically  sound  and  cost-effective  proposition.  Based  on  the  study, 
plans  were  laid  to  install  a  GPS  navigation  system  in  all  Shuttle  orbiters  beginning  in  the 
early  1980’s,  with  follow-on  goals  of  deleting  certain  equipment  from  the  baseline  navi¬ 
gation  system.  Within  about  two  years,  however,  the  decision  to  install  GPS  was 
reversed  in  favor  of  continuing  with  the  baselii»  navigation  configuration.  Certain  GPS 
provisions,  notably  antennas,  cabling,  and  bulkhead  feedthroughs,  were  nevertheless 
retained,  and  currently  exist  on  all  Orbiter  vehicles.  (Saunders,  1994,  pp.  1-13) 

The  issue  of  GPS  installation  in  the  Orbiter  fleet  surfaced  again  in  the  early 
1990’s.  Renewed  interest  was  motivated  by  the  planned  phase  out  of  TACAN  stations. 
Since  TACAN  was  used  as  the  Orbiter’s  primary  navigation  aid  following  exit  from 
blackout,  through  MSBLS  acquisition,  NASA  considered  suitable  alternatives.  Looking 
at  the  direction  the  Department  of  Defense  (DOD)  and  the  Federal  Aviation  Administra¬ 
tion  were  heading  in,  GPS  was  chosen  as  the  replacement  for  TACAN.  A  developmental 
test  for  the  Orbiter  GPS  navigation  system  flew  aboard  STS-61  in  December  1993,  and 
the  system  is  presently  expected  to  be  operational  in  1996.  (Kachmar  et.  al.,  1993,  pp. 
313-326) 
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3.  GPS  On-OTMt  DemonstradoD 

In  mid- 1992,  with  the  foundation  for  installing  an  Orbiter  GPS  navigation 
system  laid,  the  crew  of  STS-Sl  conceived  the  GOOD  DTO  as  a  low  cost  pathfinder 
project,  to  look  at  GPS  in  orbit,  to  quantify  advantages,  and  identify  potential  problem 
areas  for  Space  Shuttle  (^rations.  Data  from  the  DTO  could  then  be  used  to  comple¬ 
ment  the  more  highly  integrated  GPS  Development  Flight  Test.  Since  STS-Sl  would 
carry  another  payload  with  its  own  GPS  receiver,  the  Orbiting  Retrievable  Far  and 
Extreme  Ultraviolet  Spectrometer-Shuttle  Palette  Satellite  (ORFEUS-SPAS),  the  DTO 
would  also  permit  the  evaluation  of  relative  GPS^  Successful  utilization  of  GPS  on  the 
Orbiter  could  show  benefits  for  use  on  mher  programs,  such  as  Space  Station,  or  for  use 
as  a  utility  with  other  primary  and  secondary  payloads,  which  require  precise  location 
and  timing  information.  Initially,  goals  of  fois  experiment  were  as  follows: 

•  Evaluate  receiver  performance  in  <Ml)it  by  comparing  its  state  vector  to  that  deter¬ 
mined  by  ground  tracking  and  Orbiter  IMU’s. 

•  Evaluate  the  number  and  location  of  GPS  antennas  required  to  provide  best  naviga¬ 
tion  solutions  for  flight  deck  experiment  applications. 

•  Determine  the  quality  of  GPS  data  received  during  on-orbit  q)erations  by  collecting 
GPS  health  data. 

•  Evaluate  the  accuracy  of  relative  GPS,  using  GPS  receivers  both  in  the  crew  cabin 
and  on  ORFEUS-SPAS,  with  Orbiter  radar  and  laser  rangefinders  as  a  reference. 

•  Evaluate  postflight  the  accuracy  of  reituive  GPS  using  data  from  Orbiter  and  SPAS 
GPS  receivers. 

One  of  the  computer  displays  developed  for  this  DTO  showed  the  magnitude 
of  the  position  difference  between  the  Orbiter  GPS  and  Orbiter  IMU  based  state  vectors 
versus  time. 


^  Aspects  of  relative  GPS  are  covered  in  the  tlKsis  “Theoretical  Basis  for  State  Vector  Compari¬ 
son,  Relative  Position  Display,  and  Relative  Position/Rendezvous  Prediction”  by  LT  Lester  Makepeace, 
and  the  thesis  “NPS  State  Vector  Analysis  and  Relative  Motion  Plotting  Software  for  STS-Sl*  by  LT  Lee 
Baricer. 


4 


Since  error  in  the  IMU  based  state  vector  increased  with  time,  the  rocK  sum 
square  (RSS)  difference  (delta)  between  the  GPS  state  vector  and  IMU  state  vector  was 
expected  to  increase  with  time.  This  dilTerence  was  expected  to  collapse  to  zero  when 
Mission  Control  uplinked  a  new  state  vectm  based  on  ground  tracking.  This  behavior 
was  first  observed  in  orbit  on  flight  day  three,  when  a  20,(X)0  ft.  delta  collapsed  to  about 
3(X)  feet  following  an  update.  An  illustration  of  a  similar  event  on  the  computer  display 
is  shown  in  Figure  1  (the  x-axis  represents  time,  the  y-axis  represents  RSS  difference). 
The  close  correlation  between  expected  behavior  and  actual  behavior  indicated  die  GPS 
position  solution  to  be  near  truth. 


Figure  1:  State  Vector  Differences  vs.  Time 
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B.  GPS  OVERVIEW 

The  GPS  constellation  consists  of  21  operational  satellites,  and  three  active  spares, 
distributed  in  six  orbital  planes  with  three  or  four  operational  satellites  in  each  plane. 
The  ascending  nodes  of  each  plane  are  sq>arated  by  60°  intervals,  and  each  plane  has  an 
inclination  relative  to  the  equator  of  55°.  The  satellites  orbit  at  an  altitude  of  20,200  km, 
with  a  corresponding  period  of  12  hours.  In  comparison,  the  Space  Shuttle  orbits  at  an 
altitude  of  approximately  300  km,  with  a  corresponding  period  of  VA  hours.  The  satel¬ 
lites  are  positioned  so  that  a  minimum  of  five  will  normally  be  observable  to  a  user 
located  anywhere  on  earth. 

The  satellites  transmit  on  two  frequencies:  LI  =  1575.42  MHz  and  L2  =  1227.6 
MHz.  The  satellites  transmit  their  signals  using  spread  spectn*  a  techniques  employing 
two  different  spreading  functions:  a  1.023  MHz  coarse/acquisition  (C/A)  code  on  LI 
only  and  a  10.23  MHz  precision  (P)  code  on  both  LI  and  L2.  Both  P-code  and  C/A-code 
enable  a  receiver  to  determine  the  range  between  the  satellite  and  the  user.  Superim¬ 
posed  on  both  the  P-code  and  the  C/A-code  is  the  NAVIGATION  message  (NAV-msg), 
containing  satellite  ephemeris  data,  atmospheric  propagation  correction  data,  and  satel¬ 
lite  clock-bias  information.  The  TANS  (^uadrex  GPS  receiver  flown  on  STS-51  utilizes 
only  the  C/A-code  on  the  LI  frequency  carrier. 

Two  levels  of  navigation  are  provided  by  the  GPS;  these  are  Precise  Positioning 
Service  (PPS)  and  Standard  Positioning  Service  (SPS).  The  PPS  is  a  highly  accurate 
positioning,  velocity,  and  timing  service  which  is  made  available  only  to  authorized 
users  through  cryptographic  keys.  The  SPS  is  a  less  accurate  positioning  and  timing  ser¬ 
vice  which  is  available  to  all  GPS  users.  The  TANS  (^adrex  GPS  receiver  flown  on 
STS-51  is  an  SPS  receiver.  In  the  future,  receivers  to  be  installed  as  part  of  the  Orbiter 
GPS  navigation  system  will  be  ^PS  units.  (Kachmar,  et.al.,  1993,  pp.  313-326) 
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The  SPS  is  specified  to  (nrovide  a  100  meter  (95  %  confkieiice)  horizontal  accuracy 
to  any  GPS  user  during  peacetime.  This  is  approximately  equal  to  156  meters  three- 
dimensional  (3-D)  (95  %  confidence)  accuracy.  SPS  receivers  can  achieve  a^roximately 
337  nanosecond  (95%  confidence)  Coordinated  Universal  Time  (UTC)  time  transfer 
accuracy.  The  SPS  is  primarily  intended  for  civilian  purposes,  although  it  has  many 
peacetime  military  uses  as  well.  The  SPS  htxrizontal  accuracy  specification  include 
peacetime  degradation  of  Selective  Availability  (SA)  which  is  the  dominant  SPS  c 
source.^  The  SA  position  error  distribution  resembles  a  Gaussian  distribution  with  a 
long-term  mean  of  zero.  The  SPS  peacetime  velocity  degradation  due  to  SA  is  classified. 

The  ranging  codes  broadcast  by  the  satellites  enable  a  GPS  receiver  to  measure  the 
transit  time  of  the  signals  and  thereby  determine  the  range  between  a  satellite  and  the 
user.  The  NAV-msg  enables  a  receiver  to  calculate  the  position  of  each  satellite  at  the 
time  of  transmission  of  the  signal.  Four  satellites  are  normally  required  to  be  simulta¬ 
neously  “in  view”  of  the  receiver  for  3-D  positioning  purposes.  This  allows  the  user  3-D 
position  coordinates  and  the  user  clock  of&et  to  be  calculated  from  the  satellite  range 
and  position  data.  Treating  the  user  clock  off^t  as  an  unknown  eliminates  the  require¬ 
ment  for  users  to  be  equipped  widi  precision  clocks.  Less  than  four  satellites  can  be  used 
if  the  user  altitude  or  system  time  is  precisely  known. 

When  the  receiver  has  acquired  the  satellite  signals  from  four  GPS  satellites, 
achieved  carrier  and  code  tracking,  and  has  read  the  NAV-msg,  the  GPS  receiver  is 
ready  to  start  navigating.  The  GPS  receiver  normally  updates  its  pseudoranges  and  rela¬ 
tive  velocities  once  every  second.  The  measurements  are  termed  pseudorange  because 


^SPS  did  not  routinely  meet  accuracy  specs  while  the  GPS  system  was  undergoing  test  and  verifi¬ 
cation  by  the  DOD,  prior  to  Deconber,  19^.  SA  effects  were  often  varied  to  fiirdier  degrade  accuracy. 
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the  clock  ofbet  of  a  GPS  receiver  introduces  a  bias  to  the  true  range  of  the  satellite.  The 
GPS  receiver  must  know  the  GPS  system  time  very  accurately,  because  the  satellite  sig¬ 
nals  contain  the  time  of  transmission  from  the  satellite  in  GPS  time.  The  difference  in 
time  between  the  signal  leaving  the  satellite  and  arriving  at  the  GPS  receiver  amenna  is 
directly  proportional  to  the  range  between  the  satellite  and  the  GPS  receiver,  so  it  is  of 
the  utmost  importance  that  the  same  time  reference  is  used  by  both  the  GPS  satellites  and 
the  GPS  receive. 

The  GPS  satellites  carry  two  rubidium,  and  two  cesium  atomic  frequeiK^  stan¬ 
dards.  Howevo*,  the  GPS  receiver  is  not  required  to  have  a  hi^  accuracy  clock  such  as 
an  atomic  time  standard.  Instead,  a  crystal  oscillator  is  used  and  the  GPS  receiver  cor¬ 
rects  its  offset  from  GPS  system  time  by  making  four  pseudmrange  measurements.  The 
GPS  receiver  can  use  the  four  pseudoranges  to  solve  four  simultaneous  equations  with 
four  unknowns.  The  position  equations  are  shown  in  Figure  2.  When  the  fr>ur  equations 
are  solved,  the  GPS  receiver  has  estimates  of  its  position  and  GPS  system  time.  The 
GPS  receiver  velocity  is  calculated  using  the  same  types  of  equations,  using  relative 
velocities  instead  of  pseudoranges.  GPS  receivers  perform  most  calculations  using  an 
earth-centered  earth-fixed  coordinate  system.  They  then  convert  to  an  earth  model 
defined  by  the  World  Geodetic  System  1984  (WGS  84).  WGS  84  is  a  very  precise  model 
that  provides  a  common  grid  system  for  transformations  into  other  coordinate  systems  <n: 
map  datums. 

Satellite  coverage,  as  measured  at  the  user  antenna,  can  be  affected  by  physical 
obstructions,  vehicle  maneuvering  or  aspect,  and  basic  receiver  design.  It  therefore  can¬ 
not  be  categorized  as  a  GPS  system  requirement  or  specification.  Coverage  is  defined  by 
the  orbits  of  the  active  satellites.  The  orbits  determine  the  geometric  relationships 
between  the  satellites  and  the  user,  which  the  user  measures  as  Position  Dilution  of  Pre¬ 
cision  (PDOP).  Since  the  geometric  relationships  continuously  change  as  the  satellites 
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Rl  =  cAfj 

II 

1 

{ ;ri  -  U,  )^  +  ( yj  -  £/,  )^  +  ( Z2  -  u,  )^  =  <  «2  -  4  > 

II 

( +  ( >3  -  Uy  )^  +  ( Z3  - 1/,  =  (  Rj  -  Q, ) 

=  cAr4 

( x,  -  U,  )^  +  ( y,  -  t/y  )^  +  ( Z4  -  £/,  )^  =  ( R,  -  Ci, ) 

Rl  =  pseudorange 

c  =  speed  of  light 

=  time  difference  between  signal  leaving  the  satellite  and  arriving  at  the  receiver 


jq.  y  -v  Zi  -  satellite  position 
Uy,  =  receiver  antenna  position 
Ct,  =  receiver  clock  bias 

Figure  2:  GPS  PosithMi  Equations 

move  round  their  orbits,  so  does  the  user’s  value  of  PDOP.  PDOP  is  defined  as  the 
square  root  of  the  variances  of  the  position  errors  ( PDOP  =  ),  and 

its  effects  are  illustrated  in  Figure  3. 

C.  SUMMARY 

The  GOOD  DTO  was  designed  to  display  real  time  GPS  position  and  velocity  data 
to  the  astronaut  crew  while  on-orbit,  and  to  record  this  data  on  hard  disk  for  post-flight 
comparison  with  a  ground  generated  Best  Estimate  of  Trajectory  (BET).  It  was  also 
designed  to  display  real  time  state  vector  differences  between  the  following: 

•  TANS  GPS  and  Orbiter  IMU 

•  TANS  GPS  and  SPAS  GPS 

•  Orbiter  IMU  and  SPAS  GPS 

The  state  vector  differences  between  SPAS  GPS  and  either  the  Orbiter  IMU  or  TANS 
GPS  could  then  be  compared  against  on-board  radar  and  laser  rangefinder  results. 
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Figlire  3:  Effects  (rfnM)P 


1.  State  Vector  Difrerendiig 

Three  prerequisites  had  to  be  satisfied  before  two  state  vectors  could  be  com¬ 
pared  or  differenced.  These  were  as  fdlows: 

•  Position  and  velocity  elements  v/ect  in  die  same  frame  of  reference. 

•  Time  elements  were  in  the  same  time  scale. 

•  Time  elements  were  exacdy  matched. 

Since  the  Shuttle  navigation  system  used  the  MSO  Earth-centered  Inertial 
frame  of  reference,  and  the  GPS  receivers  utilized  the  WGS  84  Earth-centered  Earth- 
fixed  frame  of  reference,  GPS  state  vectors  were  rotated  from  the  WGS  84  reference 
frame  to  the  M50  reference  frame  prior  to  comparison  with  the  Orbiter  IMU  state  vec¬ 
tor.  This  coordinate  transformation  will  be  addressed  at  length  in  a  later  chapter. 

Since  the  Shuttle  navigation  system  used  Greenwich  Mean  Time  (GMT),  and 
the  GPS  system  used  GPS  time,  GPS  state  vector  times  were  adjusted  to  GMT  prior  to 


comparison  with  Orbiter  IMU  state  vectors.  The  Orbiter’s  clock  was  set  according  to  tiw 
National  Bureau  of  Standards  UTC  standard,  making  the  GMT  time  scale  equivalent  to 
the  UTC  time  scale  in  this  aj^lication.  The  difference  between  UTC  and  GPS  time  was 
transmitted  in  the  NAV-msg  by  GPS  satellites,  ami  this  value  was  used  for  the  adjust¬ 
ment. 

Since  the  state  vectors  being  compared  were  produced  by  independent  sys¬ 
tems,  in  general,  the  time  elements  did  not  match.  The  GOOD  DTO  utilized  Cowell’s 
method  to  “propagate”  the  earlier  state  vector  fwward  in  time  until  its  time  element 
matched  the  time  element  of  the  state  vector  it  was  being  compared  to.  Although  the  dif¬ 
ference  between  time  elements  was  typically  less  than  a  second,  it  was  significant  at 
(wbitai  velocities  of  several  kilometers  per  second.  Quick  and  accurate  propagation  of 
states  is  an  active  area  of  research,  but  will  not  be  addressed  further  in  these  pages.  ^ 

2.  State  Vector  FUtering 

Both  the  Orbiter  IMU,  and  the  SPAS  GPS  state  vector  ouq)uts  were  filtered 
in  order  to  smooth  the  ouqiuts  over  time,  and  to  improve  accuracy.  The  TANS  GPS  stat- 
evector  outputs,  however,  did  not  undergo  filtering.  Though  use  of  a  Kalman  filter  to 
smooth  the  output  would  likely  improve  the  TANS  accuracy,  a  filter  was  not  imple¬ 
mented  for  the  GOOD  DTO  due  to  time  constraints.  A  Kalman  filter  has  since  been 
designed  for  use  with  the  TANS  and  will  be  addressed  in  a  iuter  chapter. 


^See  the  master’s  diesis  of  LT  Lester  Makepeace  for  a  discussitm  of  propagation,  and  an  alternative 
to  Cowell’s  method  for  this  application. 
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n.  EXPERIMENT  DESCRlPnON 


A.  INTRODUCTION 

The  test  objectives  fo^  this  experiment  were  to  demonstrate  GI^  on-orbit  perfor¬ 
mance  at  a  relatively  low  cost.  To  meet  this  objective,  NASA  needed  to: 

•  minimize  interfaces  to  the  Orbiter; 

•  use  “off-the-sheif*  GPS  technology; 

•  designate  this  flight  test  as  a  non-cridcal  (Detailed  Test  Objective)  DTO;  and, 

•  limit  hardware/software  certification  and  qualification  to  ensuring  crew  safety. 

Other  objectives  for  this  flight  included  collecting  on-orbit  GPS  data  to  be  pro¬ 
cessed  post  flight,  demonstrating  the  GPS  performance  to  STS-Sl  crew  real  time,  and 
evaluating  potential  future  use  ftv  this  hardware  and  software. 

The  GOOD  (GPS  On-Orbit  Demonstration)  software  was  developed  for  use  on  a 
GRID  386  laptop  computer  operating  at  10  MHz.  The  desire  was  to  provide  the  crew 
with  the  capability  to  command  and  control  the  GPS  receiver,  and  to  display  and  record 
GPS  data  for  real  time  and  postflight  analysis.  In  February  1993,  when  our  Naval  Post¬ 
graduate  School  team  (LT  Lee  Barker,  LT  Les  Makepeace,  LT  Steve  Rehwald,  and  LT 
Carolyn  Tyler)  arrived  at  NASA,  Jdinson  Space  Center,  the  flight  hardware  had  already 
been  selected.  The  software  used  to  interface  with  the  TANS  GPS  receiver  was  being 
fine-tuned  to  NASA’s  needs.  Software  used  to  provide  real  time  analysis  of  the  GPS  data 
had  not  been  completed. 


B.  LEVEL  I:  HARDWARE 

The  hardware  selected  by  the  NASA  team  consisted  of  the  TANS  receiver  and  its 
associated  preamplifler,  cabling,  and  three  antennas.  Foam  spacers  were  used  to  keq) 
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the  antennas  away  from  damaging  the  sensitive  ultraviolet  radiation  protective  coating  on 
the  Orbiter’s  windows.  Radio  Frequency  (RF)  absorbers  were  placed  behiml  the  anten¬ 
nas  to  minimize  the  RF  disturbances  in  the  crew  cabin.  Finally,  the  entire  assembly  was 
attached  to  the  window  using  velcro  straps.  As  a  safety  requirement,  a  JAM  (Junction 
Adapter  Module)  was  designed  and  built.  The  JAM  was  the  only  electrical  interface  to 
the  Orbiter  and  was  required  to  protect  the  Orbiter  from  any  adverse  electrical  behavior 
generated  by  the  GPS  experiment.  The  hardware  set  up  used  aboard  Discovery  is  shown 
in  Figure  4. 


Figure  4:  Hardware  Setup 
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A  more  specific  (tescr^tion  of  the  TANS  receiver  used  in  this  e]q)mment  is  listed 


in  Table  1. 


Ihble  1:  TANS  RECEIVER  DESCRIPTION 


TANS  Description 
Code  /  Cartier  Tracked 


Channels 


Antenna  input  signals 


Position  Accuracy 


Velocity  Accuracy 


Time  Accuracy 


Dimensions: 


-  Receiver  /  Antenna 


Prime  Power 


Weight:  Receiver  /  Antenna 


Dynamic  Ct4)ability: 


-  Velocity 


-  Acceleration 


-Jerk 


Data  Interface 


Temperature: 


-  Operating  /  Non-Operating 


Altitude 


^^bration 


Shock 


Humidity 


Specification 


C/A  code,  LI 


up  to  4 


25  meters  (SEP)  without  S  A 


IQO  meters  (2DRMS)  with  SA 


0.2  m/s  without  SA 


classified  with  SA 


1  microsecond  of  UTC 


5”x9.5”x2.2”  /  3.75”x4”x0.54‘’ 


3.5  Watts  @  28  VDC 


3.5  lbs/ 0.4  lbs 


8000  m/s 


4g’s 


2g’s/sec 


RS  422  dual  channel 


-40  to  70  /  -55  to  85  degrees  C 


1100  nautical  miles 


0.04  g /Hz,  100  to  1100  Hz 


40  g  /  1 1  ms,  75  g  /  6  ms 


100%  condensing 
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The  TANS  receiver  was  chosen  because  it  was  sfnce-configured  and  it  met  the 
main  objective  of  being  available  at  a  relatively  low  cost.  A  disadvantage  to  using  this 
receiver  was  that  it  had  Standard  Positioning  Service  (SPS)  capability  and  not  the  pre¬ 
ferred  Precise  Positioning  Service  (PPS).  The  latter  c^ability  would  have  increased 
cost,  slowed  progress  and  put  special  restrictions  on  this  experiment  which  would  have 
prevented  it  from  making  STS-Sl’s  scheduled  flight  deadline.  In  the  future,  NASA 
intends  to  use  the  PPS  ct^rability  which  improves  accuracy  by  100  to  156  meters  as  com¬ 
pared  to  the  SPS  receiver.  Another  disadvantage  to  using  this  receiver  was  its  use  of  a 
deterministic  point  solution  design  instead  of  a  filter.  If  die  receiver  had  incorporated  a 
Kalman  Filter  as  a  part  of  its  design,  the  processor  unit  would  perform  calculations 
based  on  a  filtering  design  and  not  based  on  user-selected  specifications.  In  this  flight 
test,  as  an  example,  NASA  chose  a  3-D  Manual  selection,  an  option  the  TANS  provides 
the  user.  If  the  3-D  Manual  is  on,  a  three  dimensional  solution  will  not  be  calculated 
unless  four  satellites  are  in  view  and  meet  certain  requirements. 

There  were  several  settings  made  to  the  TANS  which  kept  this  receiver  from  being 
tested  in  its  best  configuration.  The  optimum  TANS  receiver/antenna  performance  was 
found  to  occur  when  antennas  were  placed  on  an  unobstructed  flat  surface,  looking 
straight  up  into  the  GPS  constellation.  Unfortunately,  the  Orbiter  does  not  fly  in  an  atti¬ 
tude  or  provide  a  window  set  up  to  facilitate  such  an  optimum  antenna  placement,  but 
rather  the  Shuttle  flies  in  a  left,  right,  or  both  wings  down  attitude  (payload  bay  down 
towards  the  earth)  with  small  restrictive  windows.  As  one  might  guess,  the  worst  posi¬ 
tion  for  the  receiver  was  when  the  Orbiter  was  in  the  latter  position,  with  all  three  anten¬ 
nas  facing  mostly  away  from  the  GPS  satellites.  Due  to  this  major  drawback,  NASA 
made  special  exceptions  to  important  settings  like  Position  Dilution  of  Precision 
(PDOP),  Signal  Level  Mask  and  Elevation  Angle. 
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The  PDOP  mask,  a  user-specified  setting,  should  normally  be  set  to  a  maximum  of 
six.  This  setting  ensures  that  four  satellites  are  in  a  satisfactory  geometrical  position. 
NASA  set  the  PDOP  at  20  because  there  was  concern  that  a  setting  of  six  might  be  too 
restrictive,  due  to  limited  viewing  angles  and  Orbiter  attitudes,  which  would  cause  the 
receiver  to  make  very  few  solutions. 

NASA  set  the  Signal  Level  Mask  specification  at  five  AMU’s  (a  linearized  mea¬ 
surement  of  noise,  usually  expressed  in  decibels)  which  is  a  very  low  setting.  Nominal 
levels  are  in  the  IS  to  25  range.  If  this  setting  had  been  lower,  more  cockpit  noise  would 
have  interfered.  As  it  was,  the  PGSC’s  (the  other  laptops  in  use)  were  picked  up  by  the 
antennas  and  tracked  by  the  receiver,  even  at  levels  below  five  AMU’s.  Trimble,  since 
this  experiment,  has  modified  the  TANS  software  package  to  correct  for  this  problem. 
The  need  for  this  low  setting  was  apparent  after  looking  at  the  way  a  satellite  was 
tracked.  In  order  to  track  the  satellite  as  long  as  possible,  the  signal  cut  off  had  to  be 
low.  Also,  the  elevation  angle  was  required  to  be  low  to  enhance  tracking,  and  in  this 
case,  NASA  set  the  angle  to  zero.  This  eiKxniraged  errOT  due  to  multipath  effects  and  the 
tendency  of  the  antenna  to  pick  up  competing  signals  (i.e.  the  PGSC)  when  tracking  the 
satellites  through  low  angles.  In  Chapter  V,  the  above  error  inducing  settings  will  be 
covered  in  greater  detail. 

C.  LEVEL  I  AND  H;  SOFTWARE  OVERVIEW 

Software  was  classified  as  either  Level  I  or  Level  II.  These  categories  reflected  not 
only  their  level  of  importance,  but  the  minimum  requirements  needed  for  the  experiment 
to  fly  aboard  the  shuttle. 

At  the  impOTtant  but  basic  first  level,  this  experiment  required  the  ability  to  oper¬ 
ate  the  TANS  receiver  and  display  certain  data  to  a  PGSC  (a  GRID  386/10)  screen.  It 
had  to  store  state  vectors  and  engineering  data  to  files,  with  the  ability  to  downlink  this 
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data.  The  downlink  capability  was  important  should  a  situation  arise  when  the  astronaut 
crew  required  troubleshooting  assistance  from  ground  experts  during  the  flight.  . 

At  the  more  advanced  level.  Level  II,  NASA  desired  some  additional  capability, 
such  as  the  ability  to  input  an  Orbiter  state  vector  either  manually  or  automatically  to 
compare  the  GPS  state  vector  measurement  to  the  Orbiter’s.  NASA  also  wanted  to  com¬ 
pare  the  TANS  GPS  measurements  with  another  GPS  receiver  used  on  the  STS-51  pay- 
load,  ORFEUS-SPAS.  One  of  the  primary  missions  for  STS-Sl  was  to  carry  this 
German-made  satellite  into  space,  release  it  to  operate  independently  for  several  days, 
and  then  rendezvous  to  recover  it  prior  to  returning  home.  The  NASA  engineers  saw  this 
as  an  opportunity  to  study  the  relative  GPS  technique. 

At  a  minimum,  NASA  wanted  to  collect  the  GPS  data  via  the  TANS  receiver  for 
postflight  studies.  Ultimately,  even  after  both  Levels  I  and  II  were  fully  developed,  this 
GOOD  test  was  only  an  experiment  (or  DTO).  It  was  to  be  operated  by  the  STS-51  crew 
on  a  not-to-interfere  basis,  only.  An  example  of  interference  occurred  during  the  rendez¬ 
vous  with  the  ORFEUS-SPAS,  which  was  one  of  the  specified  phases  of  flight  to  record 
information  for  postflight  study.  The  TANS  experiment  could  not  be  run  because  the 
antennas,  strapped  in  the  windows,  adversely  blocked  the  crew’s  view  for  rendezvous 
and  as  a  safety  of  flight  concern  interfered  with  the  crew  and  their  duties. 

D.  LEVEL  I:  SPECIFICS 

Available  in  the  Level  I  software  were  four  interface  displays.  One  display  showed 
the  user  a  current  TANS  configuration  set-up.  A  second  display  was  used  to  send  com¬ 
mands  and  requests  to  the  TANS  receiver.  A  third  showed  the  Orbiter’s  location  on  a 
world  map,  and  the  forth  display  showed  data  for  crew  monitoring.  This  data  consisted 
of  six  rows  of  information  for  six  channels  and  their  related  channel  ID,  satellite  ID, 
acquisition  flag,  ephemeris  flag,  azimuth,  elevation,  and  doppler. 


17 


Also  provided  to  the  user  was  position  and  velocity  in  two  different  coordinate  sys¬ 
tems.  The  TANS  receiver  made  measurement  calculations  using  WGS-84  .  Its  output  fcnr 
position  and  velocity  in  this  frame  was  eidier  in  cartesian  cowdinates  or  latitude,  longi¬ 
tude,  and  altitude.  A  non-trivial  coordinate  transformation  was  perfm-med  on  the  WGS- 
84  position  and  velocity  to  translate  them  into  one  of  the  key  reference  frames  used  by 
the  Orbiter,  MSO. 

E.  LEVEL  n:  SPECmCS 

The  Level  II  software  provided  real-time  graphical  displays  of  the  TANS  or 
ORFEUS-SPAS  GPS  state  vector  measurements  and  GPS/Orbiter  state  vector  compari¬ 
sons.  Figure  5  demonstrates  the  entire  GPS  configuration  in  a  block  diagram. 


Figure  5:  Experiment  Block  Diagram 


In  words,  the  TANS  GPS  information  traveled  from  the  receiver  via  the  RS-422 
connection  to  the  PGSC  715  (the  laptop  dedicated  to  GPS).  The  Orbiter  state  vector  and 
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ORFEUS-SPAS  GPS  data  traveled  from  the  PCMMU  via  another  GRID  1S3S  laptop 
running  the  PCDecom  program  through  the  RS-232  connection  to  the  PGSC  71S  to  be 
manipulated  in  the  Level  II  code.  Part  of  the  manipulation  designs  were  to  use  the  Orbit- 
er’s  GPS  (from  TANS)  and  ORFEUS-SPAS  GPS  state  vectors  in  a  rendezvous  program. 

Much  of  the  Level  II  code  was  written  by  the  NPS  team,  with  major  guidance  and 
assistance  from  a  computer  programming  wizard,  and  author  of  the  PCDecom  program, 
Mr.  Tom  Silva.  Greater  detail  about  the  flight  code  and  mathematical  derivations  is  siq>- 
plied  in  the  theses  written  by  LT  Lee  Barker  and  LT  Les  Makepeace. 
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in.  STATE  VECTOR  DIFFERENCING 

A.  DATA  ANALYSIS 

1.  Referaice  Tr^ectory 

Following  the  STS-Sl  mission,  trajectory  data  for  the  GOOD  DTO  was 
compiled  by  the  Shuttle  Navigation  group  at  NASA  Johnson  Space  Center.  In  lieu  of  a 
true  Best  Estimate  of  Trajectory,  a  reference  trajectory,  generated  by  propagating 
between  real-time  navigation  solutions  as  computed  during  the  mission,  was  created. 
The  resulting  trajectory  had  an  accuracy  of  22S  meters  and  O.IS  meters  per  second  in 
total  position  and  velocity.  This  accuracy  was  deemed  sufficient  given  the  less  than  opti¬ 
mal  test  conditions  of  the  DTO. 

Orbiter  state  vectors,  whose  time  elements  exactly  matched  those  of  GPS 
state  vectors^  were  computed  by  interpolating  within  the  reference  trajectory  using  a 
cubic  spline  technique.  This  method  introduced  an  estimated  five  meters  of  additional 
position  error  into  the  reference  data. 

2.  GPS  Data 

Most  of  the  data  output  by  the  TANS  GPS  receiver  during  the  GOOD  DTO 
was  recorded  for  postflight  analysis.  A  summary  of  recorded  data  is  shown  in  Table  2. 
Due  to  the  sheer  volume  of  data  (over  15MB)  recorded,  a  representative  sample  was 
analyzed.  Of  approximately  2500  recorded  GPS  fixes,  369  (roughly  15%)  were  chosen 
from  three  different  time  periods  for  analysis.  Selected  double  precision  XYZ,  and 
Velocity  XYZ  ECEF  binary  formatted  data  packets  were  reformatted  into  ASCII  text 
GPS  state  vectors  using  the  TANSPOST  program  supplied  by  NASA. 

^TANS  GPS  position  and  velocity  solution  times-of-fix  were  always  matched. 


TU>le2:  GPSDAIA  OUTPUT 


Packet  Description 

Packet  Size  (bytes) 

Rate 

Double  Precision  XYZ 

40 

0.33  Hz 

Double  Precision  LLA 

40 

0.33  Hz 

Velocity,  XYZ  ECEF 

24 

0.33  Hz 

Raw  measurement  output 

30-36 

0.33  Hz 

Satellite  Tracking  Status 

28 

0.33  Hz 

Satellite  Selection 

25 

0.033  Hz 

Report  Operating  Parameters 

21 

once+any  change 

Report  Control  Options 

8 

once+any  change 

Health  of  TANS 

6 

0.033  Hz 

Almanac  Information 

70 

On  request 

Almanac  Health  Page,  Hme  &  Week 

41 

On  request 

Ionospheric  Parameters 

44 

On  request 

UTC  Parameters 

43 

On  request 

Oscillator  Offset 

8 

On  request 

Satellite  Ephemeris  Status 

171 

On  request 

3.  Data  Reduction 

Before  a  comparison  could  be  made  between  the  Qrbiter  refereiK:e  state  vec¬ 
tors  and  the  TANS  GPS  state  vectors,  a  coordinate  rotation  was  required  to  transform 
the  GPS  state  vectors  from  the  WGS  84  Earth-centered  Earth-flxed  frame  of  reference  to 
the  MSO  mean  inertial  frame  of  reference.  Since  the  standard  epoch  of  B 1950.0,  upon 
which  MSO  is  based,  was  superseded  in  1976  by  the  standard  epoch  of  J2000.0,  very  lit¬ 
tle  information  appears  in  the  current  literature  for  working  with  the  MSO  system. 
Instead,  direction  is  given  for  making  transformations  to  and  from  the  J2000.0  mean 


inertial  coordinate  system.  While  NASA  continues  to  utilize  dder  transformation  theory 
to  conven  Earth-centered  Earth-fixed  coordinates  directly  to  MSO,  the  authors  chose  to 
make  the  transformation  first  into  the  J2000.0  mean  inertial  system,  and  dien  into  the 
MSO  mean  inertial  system  using  a  constant  transformation  matrix  published  by  NASA 
Jet  Propulsion  Laboratory  (Standish,  1982,  pp.  297-302)  to  relate  the  MSO  and  J2000.0 
systems.  The  transfcvmation  into  J2000.0  allowed  the  use  of  the  newer  1980  lAU  dieocy 
of  nutation  for  improved  accuracy. 

4.  Coordiiiate  Transformatioos 

The  WGS  84  system  is  mmre  accurately  referred  to  as  a  geopotential  model 
that  has  adopted  the  Conventional  Terrestrial  System  (CTS)  (1984.0),  ctefined  by  the 
Bureau  International  de  I’Heure,  as  its  reference  frame.  Earth’s  center  of  mass  is  the  ori¬ 
gin  of  the  CTS,  as  well  as  the  MSO  and  J2000.0  systems.  The  Z-axis  of  the  CTS  is 
known  as  the  (Conventional  Terrestrial  Pole  (CTP).  The  X-axis  of  die  (CTS  is  the  leto 
meridian,  and  is  used  to  derive  Universal  Time,  specifically  UTl,  in  the  same  way  the 
Greenwich  meridian  is  used  to  derive  GMT.  The  Y-axis  completes  a  right-handed  coor¬ 
dinate  system.  (DMA-TR-83S0.2,  1987,  p.  2-1) 

Transformation  of  WGS  84  coordinates  into  MSO  coordinates  requires  two 
3x3  rotation  matrices,  each  computed  for  the  time  of  the  state  vector  being  transformed. 
Transformation  of  position  vectors  actually  requires  only  a  single  rotation  matrix.  Subse¬ 
quent  discussion  will  refer  to  this  matrix  as  '‘matrix  1.”  Transformation  of  velocity  vec¬ 
tors,  however,  requires  a  second  matrix  to  account  for  the  rate  of  change  within  matrix 
1.  This  matrix  will  be  referred  to  as  “matrix  2.” 

Matrix  1  is  actually  the  product  of  five  sqiarate  3x3  matrices,  premult^lied 
to  form  a  single  matrix.  The  letters  A,  B,  C,  P,  and  M  will  be  used  to  denote  these  five 
matrices.  A  ^plies  two  rotations  for  polar  motion.  B  ^lies  a  single  rotation  for  side¬ 
real  time  (Earth  rotation).  C  applies  three  rotations  for  astronomic  nutation.  P  applies 
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three  rotations  for  general  precession^  and  M  applies  three  rotations  for  standard  epoch 
conversion.  All  matrices  with  the  exception  of  M  are  time  varying.  The  MSO  positkm 
vector  is  the  product  of  the  transpose  of  matrix  1,  and  the  WGS  84  position  vector. 

Matrix  2  is  also  the  product  of  five  sqparate  3x3  matrices.  Four  of  the  matri¬ 
ces  used  to  compute  matrix  1,  A,  C,  P,  and  are  also  used  for  matrix  2.  Unlike  the 
other  time  varying  matrices,  the  rate  of  chaise  of  the  B  matrix  is  significant;  so  another 
matrix,  denoted  by  B,  rqilaces  the  B  matrix  in  matrix  2.  B  is  the  rate  of  change  of  the  B 
matrix.  The  MSO  velocity  vector  is  computed  in  three  steps.  In  step  one,  the  product  of 
the  tranqxMe  of  matrix  2  and  the  WGS  84  position  vector  is  computed.  In  stq)  two,  the 
product  of  the  transpose  of  matrix  1  and  the  WGS  84  velocity  vecttM*  is  computed.  In 
step  three,  the  two  resulting  vectors  are  added  vecttMrially  to  form  the  MSO  velocity  vec¬ 
tor.  The  transformation  methodology  for  position  cowdinates  is  shown  in  Equation  1 
and  die  transftxmation  methodology  fm:  velocity  coordinates  is  shown  in  Equation  2. 

j  j  J 

Mm30  =  MwGS.«  '> 

X  X  X  T 

[w3m50  “  Wwgs.m 

Methods  for  computing  the  A,  B,  6,  C,  and  P  matrices  were  taken  from  The 
Astronomical  Almanac,  The  Explanatory  Supplement  to  the  Astronomical  Almanac,  and 
Defense  Mapping  Agency  Technical  Rqiort  83S0.2.  These  methods  are  summarized  in 
sub-paragraphs  a  through  d. 

a.  Polar  Motion 

Polar  motion  parameters,  Xp  and  yp,  for  the  dates  encompassing  the 
STS-Sl  mission,  were  obtained  from  the  U.S.  Naval  Observatory.  The  A  matrix  consists 
of  a  rotation  about  the  Y-axis  by  angle  -Xp^  and  a  rotation  about  the  X-axis  by  angle  -yp. 
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The  maximum  amplitude  of  these  parameten  is  approximately  0.3  arc  seconds.  This 
corresponds  tt)  about  10  meters  of  poshxHi  difference  at  Shutde  altitudes. 
b.  SidenaiTImt 

The  B  matrix  consists  of  one  rotation  about  the  Z-axis  by  an  angle  of 
A,  where  A  ==  Hq  +  AH  +oi)*(  t  -  At ).  Hq  is  Greenwich  Mean  Skfereal  Time  (GMST)  at 
0*^  UTl  on  the  day  of  interest.  Since  1984,  GMST  has  been  deflned  by  Equation  3, 
where  is  the  number  of  centuries  elapsed  since  UTl  on  2000  January  1.  The 
result  is  in  units  of  seconds  of  sidereal  time,  and  may  be  converted  to  arc  on  die  basis  of 
one  revolutkm  per  24  hours  of  sidereal  time. 

GMSTl  of  0**  UTl  =  24110.54841  +  8640184.8128667.  +  0.093 1047^  -  6.2  x  10~®7^ 

U  U  tt 

(Bq:3) 

AH  is  the  Equation  of  the  Equinoxes.  It  equals  arctan(  cose  tanAy ), 
where  e  is  true  obliquity  of  die  ecliptic,  and  Ay  is  nutation  in  longitude.  Both  e  and  Ay 
are  computed  in  the  course  of  generating  the  C  matrix.  (O*  is  the  Earth  rotation  rate  in  a 
processing  reference  fnune,  and  is  equal  to  m'  +  m,  where  co'  is  the  Earth’s  inertial  rota¬ 
tion  rate,  a  constant,  and  m  is  equal  to  7.086  x  10*^^  +  4.3  x  10'^^r„.  Time  t  is  the  time 
of  the  state  vector  being  transformed  in  seconds  since  the  beginning  of  the  day  UTC,  and 
At  is  the  difference  between  UTl  and  UTC.  Values  of  UTl  minus  UTC  were  obtained 
from  the  U.S.  Naval  Observatory  for  the  dates  encompassing  the  STS-Sl  mission.  These 
values  are  kept  below  0.7  seconds  thrcHigh  introduction  of  leap  seconds  into  UTC.  One 
second  of  time  corresponds  to  15  arc  seconds,  or  about  485  meters  of  position  difference 
at  ShutUe  altitudes. 

(1)  Change  in  Sidereal  Time  Matrix.  The  B  matrix  is  defined  as 
shown  in  Equation  4. 
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fEq:  4) 


-CD*  sinA  CO*  cosA  0 
-CO*  cosA  -CO*  sinA  0 
0  0  0 

c.  Astronomic  Nutation 

The  C  matrix  consists  of  a  rotation  about  the  X-axis  by  angle  e,  fol¬ 
lowed  by  a  rotation  about  the  Z-axis  by  ai^le  -Ay,  followed  by  a  rotation  about  the  X- 
axis  by  angle  -e.  Angle  e  is  nwan  obliquity  of  the  ecliptic,  and  is  defined  by  Equation  S, 
where  T  is  the  number  of  Julian  centuries  elapsed  since  fundamental  epoch  J2000.0  in 
barycentric  dynamical  time.  The  result  is  in  units  of  arc  seconds. 

e  =  84381.448  -  46.8157-0.000597^  +  0.0018137^  (Eq:5) 

Ay  is  nutation  in  longitude,  and  is  defined  by  Equation  6,  where  A^, 
Bi.  4ii,  a2i,  a3i,  a4i,  and  a5j  are  constants  from  the  1980  lAU  nuation  series,  shown  in 
Table  3.222. 1  of  the  Explanatory  Supplenwnt  to  the  Astronomical  Almanac,  and  (,  f,  F, 
D,  and  A  are  fundamental  arguments  of  the  1980  lAU  dieory  of  nutatioa.  The  result  of 
nutation  in  longitude  is  in  units  of  0.0001  arc  seconds. 

106 

Ay  =  5;  (Aj+BjDsinfajjl  +  ajiC  +  ajiF+a^jD  +  ajjft)  (Eq:  6) 

i  »  1 

The  fundamental  arguments  are  depicted  in  Equations  7-11.  The 
superscript  r  represents  revolutions,  and  results  are  in  units  of  arc  seconds. 

I  =  485866.733+  (1325^  +  715922.633)7+31.317^  +  0.0647^  (Eq:  7) 

r  =  1287099.804  +  (99'  +  1292581.244)  7-  0.5777^  -  0.0127^  (Eq:  8) 
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P  =  333778.877+  (1342' +295263.137)7-  13.2577^  +  0.0117^  (Eq:  9) 
D  =  1072261.307+  (1236' +  1105601.328)  7-6.8917^  + 0.0197^  (Eq:10) 
Q  =  450160.28  -  (5'  +  482890.539)  7+  7.4557^  +  0.0087^  (Eq:  1 1) 

Angle  e  is  true  obliquity  of  the  eclq>tic,  and  is  equal  to  e  +Ae,  where 
Ae  is  nutation  in  obliquity.  Ae  is  defined  by  Equation  12,  where  and  are  additional 
constants  from  the  1980  lAU  nutation  soles,  found  in  table  3.222.1  of  the  Explanatory 
Supplement  to  the  Astronomical  Almanac.  The  result  of  nutation  in  obliquity  is  in  units 
of  0.0001  arc  seconds. 

106 

Ae  =  5^  (C. +D.7)cos(a|je+a2il‘  +  a3jF+a4j7>  +  a5.a)  (Eq:12) 

i « 1 

(1)  More  Accurate  Nutation.  Av^  and  Ae^  are  corrections  to  be 
added  to  the  1980  lAU  nutations  in  longitude  and  obliquity,  Aty  and  Ae,  respectively. 
Ay^.  and  Ae,.  are  defined  by  Equations  13  and  14,  where  LQ,  OC^  and  OS^  are 
constants  from  the  corrections  to  lAU  1980  nutation  series  given  in  table  3.224. 1  of  the 
Explanatory  Supplement  to  the  Astronomical  Almanac.  The  results  of  these  nutation  cor¬ 
rections  are  in  units  of  0.(X)001  arc  seconds. 

4 

Ay^  =  2  (L5„sinA^+LC^cosA„)  (Eq:13) 

D  s  1 

4 

ACc  =  2  (ClC„cosA„  +  05„sinA^)  (Eq:14) 

n  s  1 
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An  is  equal  toa^  +  b/  +Cj^  +  +  <nQ.  where  an,  bn,  Cn,  dn, 

and  en  are  additional  constants  from  the  corrections  to  the  1980  lAU  nutation  series 
given  in  table  3.224. 1  of  the  Explanatory  Sui^lement  to  the  Astronomical  Almanac. 
d.  General  Precession 


The  P  matrix  consists  of  a  rotation  about  die  Z-axis  by  angle  -C,  fol¬ 
lowed  by  a  rotation  about  the  Y-axis  by  angle  6,  followed  by  a  rotation  about  the  Z-axis 
by  angle  -z.  Angles  C,  a^nd  z  are  defined  by  the  accumulated  precession  angles  adopted 
by  lAU  1976.  and  are  shown  in  equations  15  - 17,  respectively.  The  results  are  in  units 
of  arc  seconds. 


c  =  2306.21817+  0.301887^  +  0.0179987^ 

(Eq:15) 

z  =  2306.21817+1.094687^  +  0.0182037^ 

(Eq;16) 

0  =  2004.31097-  0.426657^-0.0418337* 

(Eq:17) 

Standard  Epoch  Conversion 

The  M  matrix  is  depicted  in  Equation  18. 

0.9999256791774783  -0.0111815116768724  -0.0048590038154553 
0.0111815116959975  0.9999374845751042  -0.0000271625775175 
0.004859003771445  -0.000027170449221  0.9999881946023742 


(Eq:18) 


5.  Summary 

The  coordinate  conversion  process  was  implemented  in  C  +  +  (Borland  ver. 
3.1).  The  program  was  run  using  an  IBM  compatible  personal  computer.  The  computer 
code  is  included  as  Appendix  A.  The  code’s  accuracy  was  validated  in  two  ways.  First, 
precession,  nutation,  and  Earth  rotation  (sidereal  time)  angles  as  shown  in  the  Astro¬ 
nomical  Almanac  were  accurately  reproduced  for  a  given  date.  Second,  a  sample  set  of 
transformed  TANS  GPS  coordinates  produced  by  the  program  was  compared  with  the 
same  set  of  coordinates,  transformed  by  the  Shuttle  Navigation  groiq)  at  NASA.  Radial 
position  differences  were  less  than  six  meters  (less  than  0.18  arc  seconds  of  rotation). 
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and  radial  velocity  differences  were  less  than  0.007  meters  per  second.  This  level  of 
accuracy  was  considered  sufficient,  given  the  differing  methods  used  in  making  the 
transformation. 

B.  NAVIGATI(»f  PERFORMANCE 

Following  coordinate  conversion,  tte  TANS  GPS  state  vectors  and  Orbiter  refer¬ 
ence  state  vectors  were  input  into  a  second  computer  program  to  compute  the  RSS  posi¬ 
tion  and  velocity  differences.  The  computer  code  utilized  is  included  as  Appendix  B. 
The  results  for  each  data  set  were  plotted,  and  the  mean  differences  and  standard  devia¬ 
tions  were  computed. 

Three  sets  of  data  were  chosen  for  analysis  of  the  TANS  GPS  receiver’s  naviga¬ 
tion  performance.  The  principle  criterion  used  in  selection  of  these  data  points  was  the 
absence  of  any  prolonged  time  interval  between  navigation  Hxes  for  the  period  under 
consideration.  For  the  periods  chosen,  time  between  fixes  is  generally  2.5  seconds,  wife 
occasional  gaps  of  up  to  7.0  seconds.  Data  recorded  during  the  STS-51  mission  was  seg¬ 
regated  into  flies  labeled  A  through  M,  O,  P,  Ql,  Q2,  Rl,  R2,  and  S  through  U.  In  gen¬ 
eral,  the  flies  corresponded  to  a  particular  event  or  activity  during  the  mission.  Large 
files  were  subdivided  into  smaller  flies  using  a  numerical  suffix,  such  as  P.001,  P.0Q2, 
etc.  The  P  flies  corresponded  to  the  crew  sleep  period  between  flight  day  6  and  7,  and 
contained  over  half  the  TANS  GPS  state  vectors  collected  during  the  mission.  All  data 
analyzed  was  taken  from  P  flies. 

Analysis  of  the  results  is  shown  in  Table  3,  and  in  Figures  6  through  11.  The  Fig¬ 
ures  depict  the  magnitude  of  the  position  and  velocity  differences,  and  show  the  range  of 
error  for  each  measurement,  as  well  as  for  the  reference  trajectory.  The  times  immedi¬ 
ately  preceeding,  and  immediately  following  the  analyzed  data  were  periods  when  the 
receiver  was  not  computing  navigation  solutions.  In  general,  a  fourth  satellite  had  just 
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TUAt  3:  COMPARISON  OF  GPS  AND  REFERENCE  TRAJECTORY 


File 

# 

of 

Samples 

Elapsed 

Time 

Position 

Difference 

Velocity 

Difference 

m:ss.s 

Mean 

(m) 

Standard 

Deviation 

Figure 

# 

Mean 

(m/s) 

Standard 

Deviation 

Figure 

# 

P.004 

98 

4:23.5 

425.8 

83.2 

6 

1.71 

0.76 

7 

P.005 

131 

5:41.0 

222.5 

30.6 

8 

0.83 

0.28 

9 

P.008 

140 

6:14.5 

322.0 

29.0 

10 

1.60 

1.02 

11 

become  usable  at  the  beginning  of  the  period  being  analyzed,  and  less  than  four  satellites 
were  usable  at  the  end  of  the  period. 

Differences  between  the  position  and  velocity  plots  for  the  different  data  sets  were 
attributable  to  a  combination  of  factors.  Among  these  were  the  number  of  usable  satel¬ 
lites  available  to  the  receiver  (4,  5,  or  6),  the  corresponding  geometry  of  the  usable  sat¬ 
ellites,  and  the  effects  of  SA.  Of  note,  SA  was  implemented  on  most,  but  not  all 
satellites  in  the  GPS  constellation,  causing  the  accuracy  of  each  measurement  to  be 
dependent  on  which  satellites  were  used  for  the  measurement. 

TANSGRAPH,  a  plotting  utility  furnished  by  NASA,  was  used  by  the  authors  to 
display  raw  data  recorded  during  the  mission.  TANSGRAPH  plots  were  generated  for 
files  P.004,  P.OQS,  and  P.008,  showing  periods  when  four  usable  GPS  satellites  were 
visible,  and  the  associated  PDOP’s  for  those  time  periods.  The  plots  are  included  as 
Appendbc  C. 

In  general,  the  position  data  was  quite  smooth,  with  the  TANS  GPS  and  Orbiter 
reference  position  error  overlapping  the  majority  of  the  time.  Conversely,  the  velocity 
data  speared  more  erratic,  with  the  TANS  GPS  and  Orbiter  reference  velocities  coin¬ 
ciding  only  a  couple  of  times  in  one  of  the  data  sets.  The  relatively  large  velocity  errors 
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would  introdiKX  significant  error  if  used  in  a  propagation  scheme,  and  would  be  unsuit¬ 
able  for  navigation  purposes  without  some  effective  form  of  filterii^.  A  few  data  points 
were  observed  to  differ  from  neighboring  points  by  an  anomalously  large  amount.  It  is 
possible  that  these  deviations  were  produced  as  an  artifact  of  data  recording,  or  data 
extraction  by  the  TANSPOST  program.  NASA  has  since  incorporated  a  six  standard 
deviatfen  rejection  scheme  into  the  GPS  processing  to  eliminate  outlying  points.  (Saun¬ 
ders,  1994,  pp.  1-13) 
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Figure  6:  File  P.004  Position  Differences 
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Figure?:  File  P.004  Vdodty  Differences 


Reference  velocity  uncertainty  (.ISm/e)  ^  Reference  velocity  uncertainty  (-.ISm/e) 
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Figure  8:  File  RODS  Position  Differences 
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Reference  position  uncertainty  C2Xm)  ^  Reference  position  uncertainty  (•230m) 
GPS  and  Reference  (P)  position  difference  ^  GPS  position  uncertainty  (-156m) 

GPS  position  uncertainty  (156m) 


Figure  9:  File  P.005  Velocity  Differences 
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Reference  velocity  uncertainty  (  ISm/s)  ^  Reference  velocity  uncertainty  (-.15m/s) 
GPS  and  Reference  (P)  velocity  difference  ^  GPS  (PPS)  velocity  uncertainty  (-.2m/8) 
GPS  (PPS)  velocity  uncertainty  (.2nV8) 


sso 


Figure  10:  File  P.008  Position  Differences 
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Reference  position  uncertainty  Reference  position  uncertainty  (•230m) 

GPS  and  Reference  IP)  position  difference  ^  GPS  position  uncertainty  (-156m) 

GPS  position  uncertainty  (156m) 


7.75 
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Reference  velocity  uncertainty  (.15m/s)  ^  Reference  velocity  uncertainty  (-.ISmfs) 

GPS  and  Reference  (P)  velocity  difference  ^  GPS  (PPS)  velocity  uncertainty  (-.2m/e) 
GPS  (PPS)  velocity  uncertainty  (.2m/e) 


IV.  STATE  VECTCMt  FILTERING 


A.  INTRODUCTION 

1.  The  Purpose  for  a  Filter 

Ch^ter  I  discussed  how  a  GPS  receiver  measures  pseudoranges  to  four  sat¬ 
ellites  in  order  to  solve  for  a  three  dimraskNial  positkm.  The  GPS  receiver  calculates  the 
user’s  position  and  GPS  time  by  knowing  the  position  of  these  four  satellites  from 
decoding  their  navigation  messages.  Pseudorange  measurements  are  made  because  the 
GPS  receiver  can  not  measure  the  exact  range  to  each  satellite.  These  measurements  ate 
ctOTupted  by  ionoqtheric  delays,  user  clock  drift,  receiver  noise  and  other  errors.  Typi¬ 
cally  a  filter  is  used  to  diaracterize  some  of  die  noise  sources  in  order  to  minimize  their 
effects  on  the  navigation  solution.  A  Kalman  Filter  can  be  used  both  widiin  the  receiver 
logic  n  in  a  post-processed  phase.  Additionally,  a  smoodiing  algmthm  could  be  implied 
in  the  post-processing  of  data.  Most  fUtering  sdiemes  studied  today  integrate  the  Kalman 
Filter  with  the  Inertial  Navigation  System  (INS)  by  using  external  informative  sources  to 
improve  position  (LORAN,  OMEGA,  laser  ranging,  etc.),  velocity  (Doppler  radar),  and 
altitude  (barometric,  radar  and  laser  altimeters).  In  this  thesis,  a  version  of  the  Kalman 
Filter  is  implemented  and  analyzed  in  the  post-processing  phase  in  which  \ua  position, 
velocity  and  GPS  time  are  known  from  STS-Sl’s  flight.  This  filter  is  an  adaptation  of  a 
Kalman  Filter  program,  written  by  Dr.  Titus,  a  professw  from  the  Naval  Postgraduate 
School,  which  uses  seven  states  including  position  (X,  Y,  Z),  velocity  (V^  ,Vy  ,V2 ), 
and  time  (t).  The  theory  of  the  Kalman  Filter  will  be  described  briefly  along  with  the 
adapted  computer  code  (see  Appendix  D).  The  final  results  will  be  analyzed  to  show  the 
effects,  both  positive  and  negative,  a  Kalman  Filtnr  has  on  this  post  flight  data. 


37 


B.  lULMAN  FILTER 


1.  Theoredcai  Model 

The  GPS  Kalman  Filter  is  a  model  of  how  its  stale  vector  is  changing  in  time 
how  the  host  vehicle  is  maneuvering  in  time.  The  state  vector  includes  parameters 
which  describe  the  moctel,  a  minimum  of  which  is  receiver  position  (X,  Y,  Z)  and  time. 
The  sinq)lest  way  to  describe  the  Kalman  Filter  is  as  a  recursive  estimatCN:  diat  {xtxhices 
a  minimum  covariance  estimate  of  the  state  vector  in  a  least  squares  sense.  The  covari¬ 
ance  estimate  or  matrix  expresses  the  statistical  uncertainty  in  the  state  vectOT.  The 
uncertainty  grows  during  long  periods  widiout  measurements.  However,  when  a  new 
measurement  becomes  available  it  will  be  weighed  heavily  regardless  of  how  noisy  it 
may  be  unless  the  filter  designer  plans  aocOTdingly. 

Figure  12,  from  the  NAVSTAR  GPS  User  Equipment  Introduction  guide¬ 
book  of  February  1991,  shows  a  simpiifmd  Kalman  Filter  diagram  illustrating  how  it 
processes  new  measurements  and  propagates  in  time. 


Figure  12:  Basic  Block  IXagnun  of  a  Kalman  Filter 
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The  basic  process  of  the  discrete  Kalman  Filter,  is  to  model  the  state  vector, 
Xt.  as  it  transitions  over  time,  from  timestq)  to  timestep.  The  matrix  is  called  the 
measurement  or  observaticm  matrix.  The  measurement  St  vectOT  is  a  fimctitMi  of  die  state 
given  by  H^.  The  following  equations  represent  two  processes.  Equation  19  for  sate  and 
Equation  20  for  measurement. 

Xt  =  0^1  •  +*,.1  (Eq.l9) 

=  Hj-Xt  +it  (Eq.20) 

The  O  matrix  is  the  transition  matrix  betv^n  the  covariance  matrix  (P|)  and  sate  (xt). 
To  account  for  the  uncertainty  in  the  state  and  measurement  models,  tla  noise  (Gaussian 
white)  terms  Ot.i  and  Xt  are  included.  The  Kalman  Filter  alternates  between  propagating 
the  state  (Xt)  and  ia  covariance  (Pj)  and  updating  these  variables  with  new  measure¬ 
ments.  The  Qt  matrix  is  the  variance  of  die  state  noise  and  it  accouna  for  the  error  in  the 
modeling  assumptions  of  O.  The  following  equations  are  used  to  express  the  propagation 
of  the  state  (x^  and  the  covariance  (Pf). 

*t  =  ^t-i  (Eq.21) 

Pt  =«l-l  Pt-l  ♦Vl  +  Ql-l  (Eq-22) 

Updating  is  defined  as  incorporating  additional  measuremena  into  the  filter 
at  regular  intervals.  The  state  immediately  iqidated  is  considered  to  be  the  most  optimal 
sate  in  the  filter.  At  this  point,  the  new  measurement  of  the  sate  is  compared  with  the 
propagated  estimate.  This  difference  is  scaled  using  the  Kalman  Gain  (Kf)  and  then  used 
in  calculating  the  new  sate  estimate.  Symbols  (-)  and  (+)  are  used  to  distinguish 
between  filter  estimates  immediately  before  and  after  a  measurement.  The  following 
equation  expresses  the  update  procedure. 

Xt  (+)  =  Xt  (-)  +  Kt'f  Zt  -  Ht  XtC-)  J 


(Eq.23) 


The  next  step  is  to  update  the  covarumce  matrix  (Pf).  In  Equation  24,  die  new  measure¬ 
ment  is  weighted  by  the  Kalman  Gain  (K|)  and  differenced  with  tte  identity  matrix  (I). 
This  result  determines  the  degree  to  which  the  covariance  matrix  is  improved  by  the  new 
nteasurement. 

Pt(+)  =  [I-KtHt]Pt(-)  (Eq.24) 

Finally,  some  understaiKling  of  the  Kalman  Gain  (K|)  is  required.  It  is  a 
result  figured  each  time  a  measurement  update  occurs.  The  calculation  is  ncM  only  based 
upon  the  propagated  covariance  matrix  of  the  previous  time,  P^  (-),  but  it  is  also  upon  the 
current  measurement  noise  covariance  (Rf)  and  the  sensitivity  of  the  measurement  to 
small  changes  in  the  state  (H^  =  SH/5xt  )* 

K,  =  P,  M-tf,  [Hi-P,  (-)  -tf,  +  R,  r‘  (Eq.25) 

To  try  to  explain  this  equation,  an  example  from  the  NAVSTAR  GPS  User 
Equipment  wUl  be  cited.  First,  assume  diat  the  state  vector  and  measurement  matrix  are 
both  in  the  same  coordinate  frame  so  that  the  Ht  matrix  becomes  the  identity  nuitrix. 
Second,  simplify  the  notation  and  matrix  formulation  to  show  the  Kalman  Gain  as  K  = 
P  /  (P  +  R),  where  P  continues  to  represent  the  covariance  and  R  represents  the  mea¬ 
surement  noise.  So  for  a  large  P,  or  uncertainty  in  the  model,  compared  to  the  uncer¬ 
tainty  in  R,  the  gain  applied  to  the  new  measurement  is  weighed  heavily  at  almost  unity. 
In  other  words,  the  propagated  state  has  too  large  of  an  uncertainty,  so  the  new  measure¬ 
ment  is  seen  as  a  better  estimate  and  is  used  as  such.  For  the  other  case,  when  there  is  a 
large  uncertainty  in  the  measurement  noise  as  compared  to  the  state  (i.e.  R>  >P),  the 
Kalman  Gain  is  very  small  and  the  new  questionable  measurement  is  weighted  by  a  small 
amount. 
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C.  AN  ADAPTATION  OF  THE  KALMAN  FILTER 
1.  Anal^rsis  of  ProbluD 

Examining  filtering  for  the  first  time,  the  decisions  concerning  what  to  filter 
and  which  parameten  to  characterize  were  the  most  challenging  aspects  to  this  problem. 
After  researchii^  Kalman  Filter  thec^  and  obtainii^  guidance  from  experts  in  this  field, 
the  problem,  at  last,  became  well  defined.  First  it  was  necessary  to  study  the  TANS 
receiver.  It  was  important  to  know  that  Trimble  had  designed  it  to  work  without  a  filter 
calculating  position  and  velocity  in  a  deterministic  manner.  The  advantage  of  accessing 
this  information  without  previous  filtering  schemes  is  that  it  allows  one  to  create  several 
filtering  designs  post-flight,  having  knowledge  of  the  data’s  behavior.  Without  access  to 
the  original  pseudorange  and  pseudorate  information,  the  problem  became  one  of  filter¬ 
ing  the  output  of  the  receiver  using  the  position  (X,  Y,  Z),  velocity  (V^,  Vy,  and 
time  (t)  parameters.  A  more  elaborate  filter  mi^t  use  other  parameters,  such  as,  GDOP, 
Carrier  to  Noise,  and  accelerations  to  assist  in  sighing  new  measurements. 

Throughout  the  entire  STS-Sl  flight,  there  were  twenty  two  periods  of 
TANS  GPS  operation.  Table  4,  on  the  following  page,  is  a  modified  table  from  “The 
First  Flight  Tests  of  GPS  on  the  Space  Shuttle.”  It  is  presented  as  a  reference  to  show 
experiment  operation  periods.  Shuttle  activity,  and  related  satellite  tracking  statistics. 

After  plotting  raw  position  and  velocity  data,  as  seen  in  Chapter  III,  it  was 
evident  that  the  position  solutions  tracked  smoothly  without  applying  filtering  tech¬ 
niques.  Choosing  one  of  the  P  files,  GPSSIP.OOS,  provided  an  acceptable  base  from 
which  to  study  the  Kalman  Filter  and  analyze  its  advantages  and  disadvantages.  How¬ 
ever,  to  carry  the  navigation  process  through  potentially  long  gaps  of  measurements 
required  a  more  elaborate  filtering  scheme  not  addressed  here. 
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Ikbie  4:  SUMMARY  OF  GOOD  DTO  PERIODS  OF  OPERATION 


Attitude 

FUe(s) 

Event/ 

Activity 

FUght 

Day(s) 

Duration 

Satellite 

Thu;king 

>3 

-ZLV,-YW 

C4).U 

- 

1,8.9 

>8hrs 

0% 

+YLV,+ZW 

EthniH 

IMU  align 

3 

>91  mins 

0  to  21.9% 

Ithru  M 

sleep 

4.5 

>8hrs 

0  to  10.7% 

+YLV,-ZW 

0 

IMU  align 

6 

0.7% 

P 

sleep 

6-7 

19.8% 

Q1  thruR2 

YLV  mnvr 

7-8 

>9hrs 

0  to  25.2% 

ACTS 

A 

ACTS 

1 

>2hrs 

0% 

Deploy 

B 

Deploy 

1 

>11  mins 

0% 

SPAS  Rdvs 

S 

SPAS 

8 

>39  mins 

0% 

T 

Rendezvous 

8 

>3hrs 

0% 

2.  Diagram  of  the  Kalman  Filter 

Figure  13  is  the  schematic  of  a  Kalman  filter  which  uses  a  signal  as  an 
input.  The  MATLAB  (Version  4.0)  code  is  included  in  Appendix  D. 


Figure  13:  Kalman  Filter  Dedgn 


42 


3.  Description  of  tiw  Katanan  FUtor  Dca^ 

The  signal  input  to  the  filter  consists  of  the  cunent  position,  velocity,  and 
time.  The  is  a  n  x  m  matrix  representing  the  Kalman  Gain.  H  is  an  m  x  n  measure¬ 
ment  or  observation  matrix  which  isolates  selected  states.  Pi^yi^  and  Pk/^.i  are  square 
matrices  [2  x  2]  representing  the  covariance  of  error  of  the  estimator  at  k  given  k  obser¬ 
vations  and  at  k  given  k-1  observations,  respectively.  The  R  matrix  (m  x  m)  is  the  cova¬ 
riance  matrix  of  the  measurement  noise  and  Q  (1  x  1)  is  the  covariance  of  the  signal 
excitation. 

The  following  formulas  are  used  in  this  Kalman  Filter  design. 

G|j  =  P|i/k-l  +  RJ*^  (Eq.26) 

^k/k  =  P  ■  G|tH]  P||/k.i  (Eq.27) 

^k+l/k  *^Pk/k^^  +  Q  (Eq.28) 

These  Kalman  Filter  equations  provide  the  gains  (G)  for  a  typical  tracking  filter  as  listed 
below. 

^^K/k  ”  ^^Wk-1  ®k  (Eq.29) 

where 

Xk/k.i  =  O  Xk.i/k.i  (Eq.30) 

H  =  [1  0] 


0  = 


1  T 
0  1 


Xk/k  = 


1  T 
0  1 


^^k-l/k-l  + 


gi(t) 

gjCt) 


(<^k  -  2^-l) 


At  the  first  observation,  when  k  =  1, 

*1/0  “  ~  82(1)  =  0 


(Eq.31) 
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At  the  second  observation,  when  k  =2,  then  gi(2)  and  g2(2)  =1/T  (noise  free). 
From  the  third  observation  on,  the  gains  will  decrease  asymptotically  to  steady  smte  val¬ 
ues  which  depend  upon  the  ratio  of  the  appropriate  term  of  the  excitation  covariance  (Q) 
and  associated  measurement  noise  variance  (R). 

D.  RESULTS  FROM  USING  A  KALMAN  FILTER 

In  the  next  few  pages,  plots  generated  using  MATLAB  are  analyzed  to  show  how 
well  this  simple  filter  performed.  Beginning  with  a  look  at  the  behavior  of  the  unfiltered 
data,  it  is  clear  that  this  series  of  data  illustrates  relatively  smooth  information.  Figure  14 
and  Figure  IS  show  the  Z  position  and  Vj  (velocity  in  the  Z  direction,  termed  Z  dot)  as 
compared  to  their  estimate,  Zhat,  generated  using  the  Kalman  Filter.  The  connecting 
line  shows  the  unfiltered  position  and  velocity  data.  The  position  plot.  Figure  14,  uses 
the  **+”  symbol  to  identify  the  filtered  position  estimate.  The  velocity  plot.  Figure  15, 
uses  the  "o”  to  discretely  show  the  trend  of  the  filter  and  one  notices  the  bias  in  this  plot 
until  Q  and  R  are  fine  tuned.  In  both  of  these  figures,  Q  =0.01  and  R  =  1 . 


RmI  vs  Predicted  in  Z  using  0^.01  and  Ral 


Figure  14:  Z  Position  •  Real  (•)  vs  Predicted  (+) 


Figure  IS:  Velocity  in  Z  Direction  •  Real  (•)  vs  Predicted  (o) 

In  the  figures  above,  one  sees  that  after  approximately  five  ii^uts,  the  filter  stead¬ 
ies  and  maintains  a  track  at  that  level.  This  characteristic  highlights  a  major  issue  con¬ 
cerning  this  filter,  initialization.  Without  initialization,  the  filter  generates  a  spike.  After 
fine  tuning  the  noise  characteristics  of  Q  and  R,  the  spike  in  some  magnitude  remains. 
Figures  16  and  17  demonstrate  the  initialization  jump  for  both  position  and  velocity.  In 
order  to  emphasize  the  spike,  these  figures  use  Q  =0.01  and  R  =1.  In  Figure  16,  the 
topmost  plot  on  the  following  page,  the  discrete  time  jumps  are  plotted  against  the  posi¬ 
tion  gain  behavior  to  show  that  the  time  steps  are  not  constant  (an  average  time  step  of 
2.5  seconds  with  seven  seconds  being  the  greatest  step)  and  correlate  the  changes  in  time 
step  with  the  gain  changes.  With  a  longer  step  than  2.5  seconds,  one  can  see  the  gain 
increase  as  it  wants  to  weigh  the  “overdue”  measurement  more  than  its  own  propagated 
state. 
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Magnitude  of  Gain 


Kalman  Gain  of  VelodtywHhQa0.01  andRsI 


Figure  17:  The  Velocity  Gain  Behavior  with  Q=0.01  and  R=1 


47 


The  unfiltered  velocity  data,  in  Chapter  III,  showed  a  mean  of  0.83  meters/second, 
which  is  unsuitable  for  navigation  purpMes.  Figure  18  shows  an  acceptable  filter  for 
position;  however.  Figure  19  reveals  a  velocity  change  which  initially  jumps  up  to  2S 
meters/second  and  then  settles  down  to  jL  5  meters/second  state,  although  stable,  it  is 
unacceptable. 


OiNefance  Plots  of  Zdot  (vaiodty}  minus  Zhatfostimat*)  using  CM).01  and  Rsi 
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figure  18:  The  Differaice  not  of  Z  minus  Zhat  (the  estimate) 


Figurel9:  The  Difference  Plot  of  Zdot  minus  Zhat  (the  estimate) 
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If  process  noise  was  zero,  Q  ^0,  and  R  was  set  arbitrarily  to  one,  the  estimaUM’ 
stability  would  converge  to  zero,  as  shown  in  Figure  19.  If  on  the  other  hand,  sensor 
noise  was  zero,  R  =0,  and  Q  was  set  arbitrarily  to  one,  the  filter  would  take  each  mea¬ 
surement,  without  weighing  them,  and  assume  they  were  correct.  Figure  20  shows  that 
there  is  no  difference  between  the  input  measurement,  Z,  and  the  estimate,  Zhat  There¬ 
fore,  the  filter  cannot  be  designed  to  remove  the  error  completely  (i.e.  Q  and  R  cannot 
equal  zero).  It  must  manage  the  disturbances,  both  sensor  and  process  noise.  The  filter’s 
accuracy  depends  upon  these  two  noise  vectors.  (Kaminar,  1993,  p.l68) 


Figure  19:  The  Position  Gain  Bdiavior  with  Q=0  and  R=1 
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RMl  v«  Prtdidtd  in  Zdol  using  0>1  snd  FW) 


Figure  21;  Vdodty  in  Z  Direction  -  Real  (-)  vs  Predicted  (o)  with  0*1  and  RsO 

The  gain,  sin^ily,  is  a  weight  emphasizing  the  incoming  measurement.  If  die  mea¬ 
surement  is  noisy,  it  should  be  de-emphasized  by  the  filter.  This  is  accomplished  by 
increasing  the  renaor  noire  <R).  However,  in  this  caw.  the  inibnnaiion  in  the  file 
appeared  to  be  steady  which  would  suggest  a  low  R.  As  for  the  process  noise  (Q).  its 
value  rqvresents  the  errors  in  the  modeling  assumptions. 

If  only  worried  about  filtering  position,  a  Q  of  0.01  proves  to  be  a  viable  setting. 
But  in  order  to  filter  the  velocity,  a  higher  Q  became  necessary.  After  running  many 
combmations  of  Q  and  R,  the  filter  began  to  provide  acceptable  velocity  variances.  The 
following  figures  illustrate  the  differences  between  conditions  of  Q  and  R  as  they  are 
changed.  Figure  22  shows  the  filtered  velociQr  differences  using  Q  =  1  and  R  =  1  which 
can  be  compared  to  Figure  23  in  which  Q  »1  and  R  =0.05,  and  finally,  the  previous 
two  plots  can  be  compared  to  the  preferred  setting  in  Figure  24,  Q  =3.5  and  R  =0.05. 
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Zhal  Zdot 


One  still  recognizes  the  spike,  albeit  its  effiKt  is  considerably  reduced  in  F^ure  24.  This 
will  only  be  corrected  with  an  initialization  tecfankpie  not  addressed  in  this  thesis. 


DMferance  Phils  of  Zdot  (valodly)  minus  ZhM(ostimate)  using  R«1  and  Qsl 


Figure  22:  The  Difllereiioe  Plot  Zdot  nsrinus  Zhat  (Qsi,  Rsi) 


Oiftarence  Plots  of  Zdot  (vslocity)  minus  Zhat(estimato)  using  R>.05  and  0^1 


Figure  23:  The  Difference  Piot  of  Zdot  minus  Zhat  (Q=l,  RsO.05) 


OllfiranM  PkMs  of  Zdot  (velocity)  minus  Zhstfsstimats)  using  fU.06  and  CM.S 


Flsure  24:  The  Difference  Plot  oi  Zdot  minus  Zhat  (Qs3^,  Rs0.05) 

E.  SUMMARY 

The  Kalman  Filter  has  many  applications  and  derivations  to  suit  ones  navigation 
requirements.  The  advantage  to  implementing  this  filter  in  software  after  the  GPS 
receiver  has  processed  the  pseudorange  and  pseudorate  information,  is  that  it  is  both  an 
inexpensive  ai^oach  to  experimental  testing  and  it  allows  access  to  raw  data,  untouched 
by  an  unproven  filter.  Without  any  previous  Shuttle  GPS  data  to  characterize,  NASA 
would  be  operating  in  unchartered  territory.  Therefore,  it  was  beneficial  to  start  at  the 
lowest  level,  the  experimental  stage,  befcnre  implementing  the  futuristic  triple  redundant 
GPS  navigation  system.  The  disadvantages  to  filtering  are  basically  the  same  as  the 
advantages  listed  above,  but  for  different  reasons.  This  relatively  inexpensive  SPS 
receiver,  if  replaced  by  a  PPS  version  with  a  filtering  scheme  using  other  characteristics 
(GDOP,  etc.)  to  weigh  incoming  measurements,  would  have  resulted  in  an  entirely  dif¬ 
ferent  thesis.  In  addition  to  the  increased  cost,  the  filter  would  have  returned  processed 
and  not  raw  data,  unfortunately,  preventing  scrnie  types  of  post-flight  analysis. 
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V.  ERR(»  SOURCES 


A.  PREDOMINANT  ERROR  SOURCES 

In  Chapters  II  and  IV,  several  errors  were  addressed  briefly;  however,  in  this 
chapter  an  exhaustive  list  will  be  discussed.  The  major  errw  in  the  navigation  solution  is 
generated  from  the  OOD’s  use  of  Selective  Availability  (SA).  System  accuracy  may  be 
degraded  by  DOD,  to  fixed  levels  by  dithering  (a  method  of  moving  the  locations  of  the 
bits  about  in  time  using  a  technique  ‘‘authorized  users”  may  remove).  The  current  policy 
implements  SA  at  a  level  giving  a  100  meter  horizontal  (two  standard  deviation)  error. 
This  intentional  error  is  the  largest  one  and  has  the  greatest  impact.  Although  other 
errors  are  addressed,  unless  a  method  such  as  differential  or  relative  GPS  is  incorpo* 
rated,  or  a  PPS  receiver  is  utilized,  these  smaller  errors  control  less  than  one  half  the 
magnitude  of  the  total  enor.  The  additional  errors  which  will  be  addressed  are  atmo¬ 
spheric  delays,  clock  differences,  q>hemeris  error,  multipath,  receiver  noise,  and  Dilu¬ 
tion  of  Precision  (DOP). 

Atmospheric  delay  is  caused  by  the  ionosphere  and  troposphere.  The  ionosphere  is 
a  layer  of  free  electrons  and  ions  above  the  atmosphere  from  approximately  100  km  to 
1000  km.  Pseudoranges  are  significantly  lengthened  because  these  charged  particles 
slow  transmissions  from  the  GPS  satellites.  The  amount  of  distortion  is  directly  prt^r- 
tional  to  the  number  of  electrons  along  the  transmission  padi.  The  error  can  range  from 
S  to  40  meters,  depending  on  the  scenario.  For  example,  a  low  elevation  angle  to  the  sat¬ 
ellite  affects  the  signal  more  than  an  overhead  view. 

Four  additional  factors  affect  the  electron  concentration  in  the  ionosphere.  They 
are  the  solar  cycle,  time  of  year,  time  of  day,  and  latitude.  The  sun  follows  an  1 1  year 
cycle  with  the  next  peak  predicted  to  occur  around  the  year  2001  to  2002  which  will 
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have  its  maximum  influence  on  the  ionosphere.  Sprite  equinox  brings  the  greatest  sea¬ 
sonal  levels  of  electron  concentration.  most  active  tinrn  of  the  day  is  at  1400  local 
time  as  compared  to  night  when  the  ionosphere  causes  its  least  intoference.  To  eliminate 
errors  due  to  the  ionosphere,  the  preferred  method  is  to  use  both  the  LI  and  L2  carrier 
frequencies,  since  its  effect  varies  with  frequency.  Use  of  both  frequencies  depends  upon 
the  receivers  capabilities.  The  TANS,  for  example,  uses  only  one.  However,  in  lieu  of 
access  to  both,  modeling  the  ionosphere  is  another  viable  option. 

The  troposphere,  represented  by  the  molecules  in  the  lower  atmosphere,  also 
iengdwns  the  transmission  path.  The  tropospheric  effect  is  proptxtional  to  the  number  of 
molecules  above  it,  which  are  separated  into  categories  of  water  and  everything  else. 
The  everydiing  else  division  is  termed  the  “dry”  category.  It  can  be  modeled  as  a  func¬ 
tion  of  elevation  angle.  The  water,  or  “wet”  ctmiponent  is  more  difficult  to  model.  The 
errors  from  these  tropospheric  effects  range  from  0.  IS  to  2.S  meters.  Their  effect  is  not 
frequency  dqpendent;  therefore,  it  is  necessary  to  use  a  model  to  decrease  the  average 
error  of  approximately  2  meters. 

Errors  caused  by  clodt  bias  and  clock  drift  seem  unavoidable.  Sources  include  the 
space  segment  (i.e.  the  GPS  satellite  clock  errors),  the  control  segment  (i.e.  the  ground 
station  clock  errors),  and  the  user  segment  (i.e.  the  receiver’s  clock  errors).  Using  four 
satellites  to  arrive  at  a  three  dimensional  solution,  the  user  clock  bias  is  determined  as 
one  of  the  four  unknowns;  therefore,  it  is  not  included  in  this  error  budget.  The  paper 
written  by  R.J.  Milliken  and  C.J.  Zoller  entitled,  “Principle  of  Operation  of  NAVSTAR 
and  System  Characteristics,”  states  that  “individual  space  vehicle  clocks,  although 
highly  stable,  may  deviate  as  much  as  976  microseconds  from  GPS  system  time.”  The 
receiver  typically  employs  the  clock  correction  coefficients  (available  in  the  navigation 
message)  in  order  to  correct  this  offset.  The  ephemeris  errors  and  the  control  segment 
clock  errors  have  a  similar  problem;  therefore,  they  can  be  addressed  together.  Briefly, 
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each  of  the  GPS  satellites  transmit  their  respective  q[>hemeridtt.  This  information  is 
updated  by  the  Master  Control  Segment  based  on  monitoring  the  individual  space  vehicle 
navigation  signals  by  four  grtnind  statioia.  This  mmhod  is  a  type  of  invorted  rai^e  pro¬ 
cess.  The  process,  albeit  a  relatively  accurate  one.  still  has  residual  uncertainties.  Thm- 
fore,  the  combined  effects  in  both  space  vehicle  clock  offsets  and  ephemeris 
determinations  are  estimated  to  be  approximately  1.5  meters.  (Milliken  and  Zol- 
ler,1980,p.  9) 

A  signal  is  not  restricted  to  following  a  direct  path,  but  may  bounce  off  electro- 
magnetically  reflective  objects.  Mult^ath  results  ffmn  having  more  than  one  propaga¬ 
tion  path  from  which  the  range  measurements  are  made.  This  error  source  is  more  likely 
to  occur  at  low  elevation  angles.  The  GPS  system  is  designed  to  minimize  the  effects  of 
multipath  by  using  L  Band  which  tends  to  diffuse,  as  a  signal,  rather  than  to  reflect.  In 
addition,  GPS  receivers  are  typically  designed  to  reject  multqmth  by  setting  signal  noise 
level  thresholds  admve  the  interfering  signals  and  by  limiting  the  elevation  angle  setting 
to  greater  than  five  degrees.  Ultimately,  multipath  error  is  estimated  to  be  1.2  meters. 

The  error  attributed  to  any  receiver  depends  upon  its  original  design,  construction, 
and  logic.  This  effect  is  termed  receiver  noise  and  resolution.  The  hardware  and  soft¬ 
ware  chosen  for  the  individual  receiver  determines  its  level  of  noise  generation.  The 
noise  can  be  generated  from  thermal  interference,  quantization  inaccuracies,  and 
dynamic  lag.  The  latter  is  due  to  the  host  vehicle’s  level  of  maneuvering.  The  estimate  of 
error  in  this  category  ranges  from  1.5  meters  to  7.5  meters. 

In  order  to  determine  the  magnitude  of  the  user  position  errors  in  the  GPS  naviga¬ 
tion  fix,  one  needs  to  combine  the  geometry  of  the  four  chosen  satellites  and  the  magni- 
mde  of  the  ranging  errors.  The  Geometric  Dilution  Of  Precision  (GDOP)  parameterizes 
the  satellite  geometry  using  PDOP  (three  dimensional  position)  and  TDOP.  For  the  pur¬ 
poses  of  this  experiment,  the  PDOP  setting  is  the  focus.  Typical  settings  for  PDOP  are 
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six  or  less.  However,  in  this  expaiinent,  this  paramo  was  allowed  to  be  an  unhealdiy 
20  in  order  to  lock  onto  four  satellites  as  often  as  possible.  It  is  difficult  to  determine 
how  much  errtx’  this  setting  induced  considering  values  greater  dian  six  are  not  written 
about  in  literature.  Using  lower  PDOPs  of  about  2.6,  produce  error  estimates  that  rai^e 
from  5.8  to  10.1  meters.  Therefore,  the  only  undisputable  conclusion  is  that  lower 
PDOPs  are  better. 

B.  ERRORS  SPECmC  TO  EXPERIMENT 

Before  this  GPS  On-Orbit  Demonstration  (GOOD)  could  be  used  by  the  STS-Sl 
flight  crew,  potential  errors  inherent  to  this  experiment  had  to  be  addressed  by  NASA 
engineers.  The  problem  areas  looked  at  and  studied  before  the  flight  were  the  effects 
from  the  Shuttle  windows,  foam  spacing  dimensions,  viewing  angle  restrictions,  cockpit 
noise  interference,  and  optimal  antenna  placement  with  respect  to  Shuttle  attitude. 

One  of  the  many  studies  performed,  prior  to  this  experiment,  looked  at  the  Shut¬ 
tle’s  window  layout  and  die  possible  effects  cause  from  the  layers  of  panes  and  the  var¬ 
ied  dimensions.  Figure  25  shows  the  locations  of  the  eight  windows,  along  with  their 
numbering  scheme.  These  window  numbem  were  used  to  identify  and  recwd  placement 
of  the  antennas  during  the  flight. 
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Each  of  the  eight  windows  has  three  panes.  These  three  panes  are  separated  by  two 
gaps  varying  in  distance  depending  on  which  window  is  addressed.  The  windows  are 
similar  in  pairs.  For  example,  windows  one  and  six  have  the  same  dimensions  and  so  on. 
The  spacing  between  the  outer  and  middle  pane  varies  because  the  crew  module  floats 
inside  the  forward  fuselage  of  the  Obiter.  Additionally,  the  size  and  shape  of  the  win¬ 
dows  vary.  The  lengths  of  the  windows  range  from  25.9  to  48.9  inches  and  the  widths 
vary  between  14.2  and  34.0  inches.  Figure  26  shows  the  dimensions  for  all  of  the  win¬ 
dows. 


Figure  26:  Window  Pane  Configuration,  Thickness  and  Size 


These  characteristics  are  important  because  in  order  to  study  the  antenna  interac¬ 
tion  with  the  windows,  the  environment  needs  to  be  fully  understood.  In  addition  to  the 
multi-pane  problem,  the  antennas  were  not  allowed  to  be  placed  directly  against  the 
Orbiter’s  window  panes  (to  avoid  damaging  the  ultraviolet  radiation  protection  coating). 
The  combination  of  both  unknowns  required  a  study  of  the  field  of  view  parameters,  sig- 


nal  strength  and  frequency  shifts.  An  actual  niock-up  was  not  available  for  studying  the 
foil  effects  that  the  window  panes  imposed  on  antenna  perfOTmance.  Instead  tht  tests 
were  simulated  as  accurately  as  possible  using  samples  of  a  thermal  pane  and  two  pres¬ 
sure  panes  together  outdoors  at  NASA’s  antenna  range.  As  Table  5  shows,  using  three 
panes  in  a  simulated  fashion,  the  best  reflected  power  percentage  and  frequency  shift 
behavior  occurred  using  three  inches  of  spacing.  The  trade  off  is  in  field  of  view.  The 
optimum  field  of  view  happens  when  the  antenna  is  pressed  against  the  window  provid¬ 
ing  a  maximum  range  of  119°  to  153°.  Due  to  die  effects  listed  below  and  knowing  the 
antenna  can  not  be  placed  at  a  zero  spacing,  the  worst  case  using  the  three  inch  spacers 
is  identified  at  a  reduced  range  of  97°  to  140°. 


Ihble  5:  EFFECTS  OF  WINDOW  GLASS  ON  ANTENNA 


Number  of 
Panes 

Spacing  between 
antenna  and  panes 
(inches) 

Reflected  Power 
(%) 

Frequency  Shift 
(MHz) 

0 

N/A 

0.01 

0 

3 

0 

79 

172 

3 

0.5 

45 

89 

3 

1 

5 

0 

3 

2 

0.04 

0 

3 

3 

0.01 

0 

The  tests  results  reported  above  provided  NASA  with  a  qualitative  assessment  of  possi¬ 
ble  behavior  which  is  attributable  to  the  simulation  restrictions  in  this  study.  As  a  final 
check,  the  GPS  antennas  were  placed  inside  Discovery  vdiile  sitting  on  the  launch  pad  at 
Kennedy  Space  Center.  The  TANS  unit  was  able  to  track  GPS  satellites  even  widi  con¬ 
siderable  blockage  of  view  due  to  the  Qrbiter’s  external  fuel  tank  and  the  launch  tower. 


Unfortunately,  position  errors  were  difficult  to  characterize  due  to  large  multipath 
effects  from  surrounding  bodies.  Although  the  three  inch  spacers  gave  higher  signal  lev¬ 
els,  STS-Sl  decided  to  use  the  one  inch  spacers  for  their  flight  in  order  to  increase  the 
field  of  view. 

An  unsuccessful  attempt  to  fly  a  similar  version  of  this  experiment  aboard  STS-56 
demonstrated  a  ‘"false-lock”  phenomenon.  This  condition  occurred  when  the  antenna/ 
preamplifier  units  were  placed  near  a  PGSC.  After  thorough  research,  NASA  discovered 
that  the  Grid  laptop  computers  were  emitting  an  interfering  signal  at  or  near  die  GPS  LI 
frequency.  It  was  also  discovered  that  this  interference  was  magnified  by  metal  enclo¬ 
sures,  such  as  the  Orbiter’s  crew  cabin.  An  RF  ^sorber  was  designed  for  this  particular 
reason,  to  counter  these  interferences  inside  the  crew  cabin.  A  maximum  size  for  these 
RF  absorbers  was  derived  from  window  viewing  requirements.  The  only  other  option 
available  to  minimize  electromagnetic  interference  was  by  setting  the  receiver’s  signal 
threshold  above  the  PGSC’s  noise  level.  The  following  figure.  Figure  27,  shows  that 
although  the  signal  threshold  of  5  AMUs  was  set,  the  TANS  receiver  continued  to  lock 
on  and  track  the  computers. 


Figure  27:  PGSC  Interference  Signal 
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Additionally,  Trimble  Navigation  modified  their  signal  processing  software;  how¬ 
ever,  it  was  unavailable  for  this  flight.  In  Figure  27,  one  can  see  the  Dof^ler  behavUx’ 
difference  between  a  “false-iock”  on  the  PGSC  and  the  true  tracking  of  a  GPS  satellite. 

C.  POST  FU(»T  ANALYSIS 

The  Shuttle  operated  in  three  different  attitudes  during  this  experiment.  These  are 
shown  in  Figure  28.  Of  the  three  attitudes,  the  worst  condition  was  Attitude  3.  The 
antennas,  all  feeing  toward  the  earth,  were  never  able  to  get  four  satellites  in  view  in 
order  to  calculate  a  three  dimensional  solution.  In  fact,  the  best  attitude  of  the  three 
tested  was  Attitude  1.  The  “P”  series  of  files  were  obtained  during  a  sleep  period  in  Atti¬ 
tude  1  between  flight  day  six  and  seven  and  were  used  extensively  in  post  flight  analysis. 
For  this  period,  the  three  antennas  were  placed  in  windows  seven,  six,  and  four.  Out  of 
approximately  2500  samples,  1464  were  found  in  the  “P”  series  of  files. 


Figure  28:  Orbiter  Attitudes  During  GOOD  DTO  Operations 
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In  conclusion,  the  attitude  of  the  Orbiter  provided  a  suboptimal  environnwnt  for 
viewing  the  GPS  constellation  given  that  the  antennas  were  required  to  be  inside  the 
crew  cabin.  Placing  the  antennas  outside  in  the  Orbiter’s  Payload  Bay,  for  example, 
would  improve  the  performance  of  the  receiver  because  the  issues  concerning  foam 
spacing  and  cockpit  noise  would  disappear.  However,  the  field  of  view  and  other  poten¬ 
tial  errors  would  need  to  be  characterized.  Finally,  the  errcHS  discussed  in  this  chapter 
were  recognized  by  NASA  and  future  follow-on  experiments  will  be  modified  to  mini¬ 
mize  these  problems  areas. 
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VI.  CONCLUSIONS 


A.  FLIGHT  HARDWARE 

One  of  the  chief  aims  of  the  GOOD  DTO  was  to  keq>  costs  to  a  minimum.  This 
goal  was  readily  achieved,  but  at  the  expense  of  some  capabilities  that  might  otherwise 
have  been  available.  Two  specific  improvements  in  flight  hardware  would  greatly 
enhance  the  performance  of  GPS  aboard  the  Space  Shuttle.  Each  of  these  improvements 
would,  however,  have  a  significant  cost  associated  widi  it. 

The  first  suggested  improvement  is  the  use  of  space  qualified  equipment.  This  fea¬ 
ture  would  permit  components  to  be  located  outside  die  fnendly  environment  of  the  crew 
cabin,  in  the  Shutde’s  payload  bay.  This  is  particularly  desiridile  in  the  case  of  die  anten¬ 
nas,  whose  in-cabin  field  of  view  is  gready  restricted.  Bettor  visibility  should  enable  four 
satellites  to  be  in  view  for  a  much  higher  percentage  of  time,  resulting  in  more  consistent 
availability  of  the  GPS  state.  A  corresponding  reduction  in  PDOP  would  also  be 
expected,  correlating  to  improved  GPS  navigation  accuracy. 

The  second  suggested  improvement  is  the  use  of  a  Precise  Positioning  Service 
receiver.  This  is  particularly  desired  because  of  the  unsatisfactory  velocity  results 
obtained  by  this  experiment  due  to  Selective  Availability  (SA) .  The  velocity  accuracy  of 
the  PPS  is  0.2  meters  per  second.  This  is  the  order  of  magnitude  of  the  desired  accuracy 
for  velocity  information.  The  velocity  accuracy  of  the  Standard  Positioning  Service  is 
classified,  but  errors  observed  during  this  experiment  were  on  the  order  of  2  meters  per 
second.  Velocity  errors  this  large  are  unsuitable  for  use  in  navigation,  as  even  small 
velocity  errors  cause  large  position  errors  when  used  for  propagating  states.  Use  of  a 
Kalman  Filter  earlier  in  the  receiver’s  logic  would  smooth  both  position  and  velocity  to 
an  acceptable  level.  While  a  filtering  scheme  may  succeed  in  reducing  some  of  the  error. 
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the  only  way  to  consistently  achieve  desired  velocity  accuracies  is  through  use  of  a  PPS 
receiver. 

B.  FUGHT  SOFTWARE 

The  flight  software  functioned  as  it  was  designed,  successfully  recording  data  for 
postflight  analysis.  The  only  desirable  modific^on  to  the  software  would  be  the  addition 
of  a  relatively  fast  filtering  algorithm,  such  as  a  version  of  the  Kalman  Filter  designed  to 
smooth  GPS  state  vector  output.  Again,  incorporation  of  this  feature  would  carry  an 
associated  cost  of  further  integration  and  testing. 

C.  FUTURE  APPUCATIONS 

Each  of  the  suggested  improvements  to  the  GOOD  DTO  are  already  incorporated 
as  part  of  the  design  of  the  GPS  navigation  system  being  designed  for  the  Shuttle. 
(Kachmar,  1993,  pp.  313-326).  No  plans  exist  to  fly  the  GOOD  DTO,  as  configured  on 
STS-51,  again  at  a  future  date.  However,  portions  of  the  experiment  have  been  utilized 
f(x  two  Shuttle  missions  since  the  September  1993  STS-Sl  flight,  and  STS-66  is  sched¬ 
uled  to  carry  a  portable  GPS  receiver,  with  antennas  mounted  in  the  payload  bay,  later 
this  year.  The  GOOD  DTO  has  shown  that  good  on  orbit  position  accuracy  can  be 
obtained  using  an  inexpensive  portable  GPS  receiver.  The  output  of  such  a  receiver  can 
be  used  by  any  application  or  experiment  requiring  precise  timing  or  position  informa¬ 
tion. 
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AFFENDK  A 


/*  */ 

/*  This  program  trauisforms  a  state  vector  from  the  ECEF_WGS_84  to  the*/ 

/*  ECI_M50  cartesian  reference  frame.  */ 

/*  */ 

/*  Author:  LT  Stephen  P.  Rehwald,  USN  */ 

/*  Date:  06  March  1994  */ 

/*  */ 

/*  Functions  from  "Numerical  Recipes  in  C",  Press,  W.H.,  et  al,  */ 

/*  Cambridge  University  Press,  1988  were  utilized  in  this  program.  */ 

/*  */ 

/*  (As  written,  this  code  is  valid  for  state  vectors  having  time  */ 

/*  elements  between  MJD  49247.0  and  MJD  49248.0,  as  this  fulfilled  the  */ 
/*  author's  needs.  Other  times  may  be  used  by  incorporating  */ 

/*  corresponding  lERS  polar  motion,  UT1_UTC__0FPSET,  and  applicable  leap  */ 
/'*  second  adjustments.  Slight  modification  of  the  routines  for  */ 

/*  confuting  UTl  time  and  for  interpolating  polar  motion  data  will  be  */ 
/*  required  to  reflect  the  newly  incorporated  data.)  */ 

/*  */ 


/  It  ******  If*******************  *********************  It*  ******  ******  ****11  I 

#define  PI  3.1415926535897932385 
#define  DEGjrO_RAD  (PI/180.0) 

#define  ARCSEC_TO_RAD  (PI/ (180 . 0*3600 . 0) ) 

#define  SEC_T0_RAD  ( (2 . 0*PI) /8640r, .  0) 

#define  ONE_REV  1296000.0 

#define  DAY  86400.0 

#define  CENTURY  36525.0 

#define  OMEGA_PRIME  .000072921151467 

!************************************************************************! 
/*  NUMBER_OF_VECTORS  is  the  number  of  position/velocity  vectors  to  */ 

/*  process  from  the  files  used  as  input  to  this  program.  Each  file  */ 

/*  must,  at  a  minimum,  contain  this  number  of  vectors.  */ 

/♦♦***♦**********★******♦*♦*************♦************♦********♦**********/ 

#define  NUMBER_OF_VECTORS  139 

/♦*****************************************  *************************  j 

/*  x_P_49247  and  Y_P_49247  are  the  euigular  c.  ^acements  in  arc  seconds  */ 
/*  of  the  Celestial  Ephemeris  Pole  (CEP)  from  the  Conventional  */ 

/*  Terrestial  Pole  (CTP)  (CTP=ECEF_WGS_84  Z  socis)  in  effect  at  Modified  */ 
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/*  Julian  Date  (MJD)  49247.0.  Likewise,  X_P__49248  and  Y_P_49248  are  */ 
/*  these  values  at  MJD  49248.0.  UT1_UTC_49247  and  UT1_UTC_49248  are  the  */ 
/*  differences  in  seconds  between  the  UTl  and  UTC  timescales  at  MJD  */ 
/*  49247.0  and  MJD  49248.0  respectively.  Ex:  UT1-UT1_UTC_49247=UTC.  */ 
/*  TAI_UTC_OFFSET  is  the  integer  difference  in  seconds  between  */ 
/*  International  Atomic  Time  (TAD  and  UTC.  Ex:  TAI-TAI_UTC_OFFSET=UTC.  */ 
I*  This  difference  is  exact,  and  changes  at  not  less  than  six  month  *! 
/*  intervals  through  introduction  of  leap  seconds.  ‘IT)T_TAI_OFFSET  is  */ 
I*  the  set  difference  in  seconds  between  Terrestrial  Dynamical  Time  */ 
/*  (TDT)  and  TAI.  Ex:  TDT-TDT_TAI_OFFSET»TAI .  Future  predicted,  cuxd  */ 
/*  past  observed  values  of  x_p,  y_p,  &  utl_utc,  as  well  as  */ 
/*  TAI_UTC_OFFSET  and  TDT_TAI_OFFSET  are  published  by  the  International  */ 
/*  Earth  Rotation  Service  (lERS)  (in  the  U.S.,  by  the  National  Earth  */ 
/*  Orientation  Service  (NEOS) ) .  The  following  were  obtained  by  */ 
/*  directing  an  cuionymous  FTP  into  MAIA.USNO.NAVY.MIL  to  access  the  */ 
/*  SER7  directory.  */ 


/************************************************************************/ 

#define  X_P_49247  -.097 
#define  Y_P_49247  .345 
#define  UT1_UTC_49247  .4570 
#define  X_P_49248  -.097 
#define  Y_P_49248  .347 
ttdefine  UT1_UTC_49248  .4543 
#define  TAI_UTC_OFFSET  28.0 
#define  TDT_TAI_OFFSET  32.184 

/************************************************************************/ 
/*  GPS  UTC  OFFSET  is  the  inte'  difference  in  seconds  between  GPS  time  */ 


/*  and  Coordinated  Universal  (UTC) .  Ex:  GPS-GPS_UTC_OFFSET=UTC.  */ 
/♦  For  Standard  Positioning  Service  (SPS) ,  GPS  time  is  accurate  to  */ 
/*  within  337  ncuioseconds  of  UTC  time  after  GPS_UTC_OFFSET  is  applied.  */ 
I*  GPS_WEEK_NUMBER  is  the  number  of  weeks  elapsed  since  06  JAN  80. 

/*  GPS_UTC_OFFSET  and  GPS_WEEK_NUMBER  are  transmitted  in  NAV-messages  */ 
/*  from  GPS  Satellite  Vehicles  (SVs)  .  JD__GPS_WEEK_0  is  the  Juliaui  Date  */ 
/*  at  Oh  06  JAN  80  in  UTC.  */ 


/**♦***♦****♦********♦******************♦♦******♦♦********************♦♦*/ 

#define  GPS_UTC_OFFSET  9.0 
#define  GPS_WEEK_NUMBER  714.0 
#define  JD_GPS_WEEK_0  2444244.5 

/*  JD_00_JAN_93  is  the  Julian  Date  at  Oh  00  JAN  93  in  UTC.  JD_J2000  is  */ 
/♦  the  Julian  Date  at  Standard  Epoch  J2000.0  in  Barycentric  Dynamical  */ 
/♦  Time  (TDB) .  */ 
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/ 


/ 


#define  JD_00_JAN_93  2448987.5 
#de£ine  JD_J2000  2451545.0 

#include  <bcd.h> 

#include  <inath.h> 

#include  <stdio.h> 

#include  <stdlib.h> 

void  nerror (char*) ; 

FILE  *j,*k; 
int  n,o,i; 
long  double 

a_i[106]  [5]  ,s_iI1061  [2]  .c_i(106]  [2]  .a_nl4l  [5]  ,s_n(4]  [2]  ,c_n(4l  [2]  ; 

long  doiible  m_0[3I  [31 ,  t_0  (NUMBBR_OP_VBCTORSl  [6]  ; 

long  double  v_0  [NOMBER__OF_VECTORS]  [61  ,  v_l  [NUMBER_OF_VECTORSl  [61  ; 

V 

V 

*/ 

void  nerror (char  error_text [ 1 ) 

{ 

fprintf  (stderr,  "Ntiraerical  Recipes  run-time  error. .. \n")  ; 

f print f (stderr, "%8\n" , error_text) ; 

fprintf (stderr, " . . .now  exiting  to  system. . .Xn") ; 

_exit ( 1 ) ; 

} 


/* 

/*  nerror  is  Numerical  Recipes  standard  error  hamdler 
/* 


/* ♦******♦**♦**♦***********♦******♦♦***♦♦*******♦*********♦********♦*****/ 
/* 

/*  a_i[]  [l,s_i[l  [1,  and  c_i[l  [1  are  the  lAU  1980  Nutation  series 
/*  multipliers  and  coefficients  (Table  3.222.1,  Explanatory  Supplement 
/*  to  the  Astronomical  Almanac,  1992) . 

/*  a_n[]  [l,s_n[l  [1,  cuid  c_n[l  [1  are  multipliers  cuid  coefficients  for 
/*  corrections  to  the  lAU  1980  Nutation  series  (Table  3.224.1, 

Explcuiatory  Supplement  to  the  Astronomical  Alamanac,  1992) . 
m_0[l  [1  is  the  treuisformation  matrix  from  B1950.0  to  J2000.0  a.k.a. 
DE118/LE62  to  DE200/LE200  (E.M.  Standish,  Astronomy  and  Astrophysics 
114,  pp.  297-302,  1982). 

t_0[il  [01  is  the  time  element  of  the  ith  state  vector  in  seconds 
/♦  since  the  beginning  of  the  week  (GPS  time)  .  t__0[il  [01«t_0[il  [11  . 

/*  t_0[il  [21  is  the  time  element  of  the  ith  state  vector  in  hours 
/*  since  the  beginning  of  the  year  (UTC)  (Consistent  with  the  Orbiter 


/* 

/* 

/* 

/* 

/* 


*/ 

*/ 

*/ 

*/ 

*/ 

*/ 

*/ 

*/ 

V 

*/ 

*/ 

*/ 

*/ 

*/ 
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/*  reference  trajectory  state  vectors  supplied  by  NASA  JSC) .  */ 

/*  t_0[i]  [3]  is  the  Julian  Date  of  the  ith  state  vector  in  TDB.  */ 

/*  t__0[i]  [4]  is  the  Julian  Date  at  Oh  on  the  day  of  interest  in  un.  */ 

f*  t_0[i]  [5]  is  the  time  el«iient  of  the  ith  state  vector  in  seconds  */ 

/*  since  the  beginning  of  the  day  (UTC) .  */ 

/*  v_0ti]  [0],  v_0[i]  [1],  and  v_0[il  [2]  are  the  X,  Y,  and  Z  cartesiam  */ 

f*  position  coordinates,  respectively,  of  the  ith  state  vector  in  the  */ 

f*  ECEF_WGS_84  reference  frame.  */ 

/*  v_0[i]  [3],  v_0ti]  [4],  and  v__0[i]  [5]  are  the  X_DOT,  Y_DOT,  and  Z_DOT  */ 
f*  cartesiam  velocity  coordinates,  respectively,  of  the  ith  state  */ 

/*  vector  in  the  ECEP_WGS_84  reference  frame.  */ 

f*  v_l[il  [01,  v__l[i]  [1],  auid  v_l[i]  [2]  are  the  X,  Y,  auid  Z  cartesian  */ 

f*  position  coordinates,  respectively,  of  the  ith  state  vector  in  the  */ 

/*  ECI_M50  reference  frame.  */ 

/♦  v_lti]  [2],  v_l[i]  [4],  and  v_l[i]  [5]  are  the  X_DOT,  Y_DOT,  auid  Z_DOT  */ 
/*  cartesian  velocity  coordinates,  respectively,  of  the  ith  state  */ 

f*  vector  in  the  ECI_M50  reference  frame.  */ 

t*  */ 

/*  ECEF_W6S_84  coordinates /GPS  time  are  read  from  ASCII  position  and  */ 

(*  velocity  data  files,  created  by  the  TANSPOST  prograun,  which  extracts  */ 

f*  desired  data  from  the  binary  data  files  created  by  the  TANSIO  */ 

f*  program  on  STS-51.  Transformed  ECIJMSO  coordinates /DTC  are  written  */ 

f*  to  an  ASCII  data  file.  */ 

(*  FILE  SPECIFICATIONS:  -ECEF_WGS_84  Position  File  Naune  =  statep.OOl  */ 
(*  -BCEF_WGS_84  Velocity  File  Naune  =  statev.OOl  */ 

I*  -ECI_M50  State  Vector  File  Name  =  stateout.OOl  */ 

f*  -Data  appearing  on  the  same  line  in  both  input  */ 

f*  files  must  correspond  to  the  same  time  */ 

f*  element.  */ 

f*  -The  number  of  vectors  in  the  input  files  must  */ 

I*  be  equal  to  or  greater  thaui  the  number  */ 

f*  appearing  in  the  #define  NUMBER_OF_VECTORS  */ 

f*  preprocessor  statement.  */ 

I*  */ 

(*  From  TANS  packet  41:  */ 

f*  GPS_UTC_OFFSET  =  9  seconds  throughout  the  mission.  */ 

f*  GPS_WEEK_NUMBER  »  714  Starting  12  SEP  94  */ 

f*  GPS_WEEK_NUMBER  =  715  Starting  19  SEP  94  */ 

'*  The  correct  numbers  must  appear  in  the  #define  preprocessor  */ 

statements.  */ 

n,  */ 

'*  Coordinate  conversion  methodology:  */ 

ECI_M50  position  =  [ed)cpm_0]  trauispose  *  ECEF_WGS_84  position  */ 

'*  ECI_M50  velocity  =  [adD_dotcpm_0]  transpose  *  ECEF_WGS_84  position  +  */ 

'*  [abcpm_0]  trauispose  *  ECEF_WGS_84  velocity  */ 

'*  a,  b,  b_dot,  c,  p,  auid  m_0  are  rotation  matrices  for  polar  motion,  */ 
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/*  Barth  rotation,  rate  of  change  of  Barth  rotation  matrix,  astronomic  */ 
I*  nutation,  general  precession,  and  steuidard  epoch  conversion,  */ 
/*  respectively,  a  through  p  are  time  varying  3x3  matrices,  thoigh  only  */ 
/*  the  rate  of  change  of  b  is  significant.  m__0  is  a  constant  3x3  */ 
/*  matrix.  */ 

/*  V 


void  main (void) 

{ 

long  doiible  longestsO.O; 

//  Initialize  Nutation  series  a_i[]  [1  s_i[]  [1  c_i[]  [1  a__n[]  []  s_n[]  [1  c_n(]  []* 
for (na0;n<106;n++) 

{ 

for (O*0;O<5;O++) 

{ 

a_i [n] [o] =0.0; 

} 

for (O»0;O<2/O++) 

{ 

s_i (nj  (oj  =0.0; 
c_i[n]  [ol=0.0; 

} 

} 

for (n=0 ; n<4 ; n++ ) 

{ 

for (o=0;O<5;o++) 

{ 

a_n [n] (o] =0.0; 

} 

for ( 0=0 ; 0<2 ; 0++ ) 

{ 

s_n[n] [o] =0.0; 
c_n [n] [o] =0.0; 

} 

} 

a_i[0]  [4]  =1.0;  a_i[l]  [4]  =2.0;  a_i[2]  [0]=-2.0;  a_i[2]  [21=2.0;  a_i[2]  [4]=1.0; 

a_i[3][0]=2.0;  a_i[3]  [2]=-2.0;  a_i  [4]  [0]  =-2 . 0 ;  a_i[4][2]=2.0; 

a_i[4][43=2.0;  a_i[5]  [01=1.0;  a_i [51  [11  =-l . 0 ;  a_i[5l[31=-1.0; 

a_i[61  [ll=-2.0;  a_i  [61  [21=2.0;  a_i  [61  [31  =-2 . 0 ;  a_i[61[41=1.0; 

a_i  [71  [01=2.0;  a_i[71  [21=-2.0;  a_i  [7]  [41  =1 . 0 ;  a_i[81[21=2.0; 

a_i[81  [3l=-2.0;  a_i  [8l  [4l  =2 . 0 ;  a_i  [91  [11  =1 . 0 ;  a_i[l01[ll=1.0; 

a_i  [101  [21=2.0;  a_i[101  [31=-2.0;  a_i  [101  [43  =2 . 0 ;  a_i[lll[ll=-1.0; 

a_i  [111  [21=2.0;  a_i[lll  [31=-2.0;  a_i  [111  [41  =2 . 0 ;  a_i[121[21=2.0; 

a  i[12l  [31=-2.0;  a_i [121  [41 =1 . 0 ;  a_i [131  [01 =2 . 0 ;  a_i[131[31=-2.0; 

a^i [141  [21=2.0;  a_i [141  [31 =-2 .0;  a  i [151  [11 =2 . 0 ;  a_i [161  [11=1.0; 
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a_i[16]  [41-1.0;  a_i[171  [11-2.0; 
a_i[171[4]»2.0;  a_i[181  [11— 1.0 
a_i[191  [31-2.0;  a_i[191  [41-1.0; 
a_i[201  [31— 2.0;  a_i  [201  [41-1.0 
a_i[211  [41=1.0;  a_i[221  [11-1.0; 
a_i[221  [41=1.0;  a_i[231  [01=1.0; 
a_i[24l  [11=1.0;  a_i[24l  [31=-2.0 
a_i[25l  [41=1.0;  a_i[261  [11=1.0; 
a_i[27l  [11=1.0;  a_i[271  [41=2.0; 
a_i[281  [41=1.0;  a_i[291  [11=1.0; 
a_i[301  [21=2.0;  a_i[301  [41=2.0; 
a_i[321  [41=1.0;  a_i[331  [01-1.0; 
a_i[34]  [01  =1.0;  a_i[34l  [31--2.0 
a_i[351  [41=2.0;  a_i[361  [31=2.0; 
a_i[381  [01  =-1.0;  a_i  [381  [41=1.0 
a_i[391  [31=2.0;  a_i[391  [41=2.0; 
a_i[40l  [41=1.0;  a_i[411  [21=2.0; 
a_i[421  [01=2.0;  a_i[431  [01=1.0; 
a_i[43l  [41=2.0;  a_i  [441  [01=2.0; 
a_i[45l  [21=2.0;  a_i[461  [01--1.0 
a_i[471  [01=-1.0;  a_i[47l  [31-2.0, 
a_i[48l  [31—2.0;  a_i[48l  [41=1.0, 
a_i[491  [31=2.0;  a_i[491  [4J=1.0; 
a_i[501  [31—2.0;  a_i[511  [11  =1.0 i 
a_i[521  [11  =-1.0;  a_i[521  [21=2.0; 
a_i[531  [21=2.0;  a_i[531  [31=2.0; 
a_i[54l  [31=2.0;  a_i[551  [01=2.0; 
a_i[55l  [41=2.0;  a_i[561  [31-2.0; 
a_i[57l  [31=2.0;  a_i[571  [41=1.0; 
a_i[581  [31  =-2.0;  a_i  [581  [41=1.0; 
a_i[60l  [01=1.0;  a_i[601  [11  =-1.0; 
a_i[611  [41=1.0;  a_i[621  [11=1.0; 
a_i[63l  [21--2.0;  a_i[64l  [31=1.0; 
a_i[661  [01=1.0;  a_i[661  [21=2.0; 
a_i[67l  [21=2.0;  a_i[671  [41=2.0; 
a_i[681  [21=2.0;  a_i[681  [31=2.0; 
a_i[691  [41=1.0;  a_i[701  [01=3.0; 
a_i[711  [11  =-1.0;  a_i[711  [21=2.0; 
a_i[72l  [01=1.0;  a_i[721  [11=1.0; 
a_i[731  [01  =-1.0;  a_i[731  [21=2.0; 
a_i[74l  [01=2.0;  a_i[741  [41=1.0; 
a_i[761  [01=3.0;  a_i[771  [21=2.0; 
a_i[781  [01 =-1.0;  a_i[781  [41=2.0; 
a_i[80]  [01  =-2.0;  a_i[801  [21=2.0; 
a_i[811  [01--1.0;  a_i[811  [21=2.0; 
a_i[821 [01=2.0;  a_i[821 [31 =-4.0; 


a_i[171  (21-2.0;  a_i[17l  [31—2.0; 

;  a_i [181 (41-1.0;  a_i [191 [01 =-2 . 0 ; 
a_i[201  [11—1.0;  a_i[201  [21=2.0; 

;  a_i  [211  [01-2.0;  a_i  [211  [31  =-2 . 0  ; 
a_i[221  [21=2.0;  a_i[221  [31  =-2.0; 
a_i[231  [31  =-1.0;  a_i[24l  [01=2.0; 

;  a_i[251  [21=-2.0;  a_i [251  [31 =2 . 0 ; 
a_i[261  (21—2.0;  a_i[261  [31=2.0; 
a_i[281  [01—1.0;  a_i[281  [31=1.0; 
a_i[291  [21-2.0;  a_i[291  [31—2.0; 
a_i[311  [01=1.0;  a_i[321  [21=2.0; 
a_i[331  [21=2.0;  a_i[331  [41=2.0; 

;  a_i[351  [01  — 1.0;  a_i [351  [21 =2 . 0 ; 
a_i[37l  [01=1.0;  a_i[37l  [41=1.0; 

?  a_i[391  [01  =-1.0;  a_i[391  [21=2.0; 
a_i[401  [01=1.0;  a_i[401  [21=2.0; 
a_i[411  [31=2.0;  a_i  [411  [41=2.0; 
a_i[431  [21=2.0;  a_i[431  [31  =-2.0; 
a_i[44l  [21=2.0;  a_i[441  [41=2.0; 
a_i[461  [21=2.0;  a_i  [461  [41=1.0; 
a_i[471  [41=1.0;  a_i[481  [01=1.0; 
a_i[491  [01  =-1.0;  a_i[491  [21-2.0; 
a_i(50J  (OJ-l.O;  a_i  [501  [11=1.0; 

'  a_i  [511  [21=2.0;  a_i  [511  [41  =2 . 0  ; 

■  a_i [521 [41=2.0;  a_i [531 [01 =1 . 0 ; 
a_i[531  [41-2.0;  a_i[54l  [01=1.0; 
a_i[551  [21=2.0;  a_i[55l  [31  =-2.0; 
a_i[561  [41=1.0;  a_i[57l  [21=2.0; 
a_i[581  [01=1.0;  a_i[581  [21=2.0; 
a_i[591  [31  =-2.0;  a_i[591  [41=1.0; 
a_i  [611  [01=2.0;  a_i[611  [21  =2.0; 
a_i[621  [31=-2.0;  a_i[631  [01=1.0; 

a_i[651  [01=1.0;  a_i[65l  [11=1.0; 
a_i[671  [01=1.0;  a_i[67l  [11  =-1.0; 
a_i[681  [01 =-1.0;  a_i[681  [11 =-1.0; 
a_i[681  [41=2.0;  a_i[691  [01  =-2.0; 
a_i[701  [21=2.0;  a_i[70l  [41=2.0; 

a_i[711  [31=2.0;  a_i[711  [41=2.0; 
a_i[721  [21=2.0;  a_i[721  [41=2.0; 

a_i[731  [3]  =-2.0;  a_i[731  [41=1.0; 
a_i[751  [01=1.0;  a_i[75l  [41=2.0; 
a_i[771  [31=1.0;  a_i[77l  [41=2.0; 
a_i[791  [01=1.0;  a_i[791  [31  =-4.0; 
a_i[801  [31=2.0;  a_i[801  [41=2.0; 
a_i[8ll  [31=4.0;  a_i[811  [41=2.0; 
a_i[831  [01=1.0;  a_i[831  [11=1.0; 
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a_i [83]  [21-2.0;  a_i [83]  [3] — 2 .0 ;  a_i[83]  [4] -2.0;  a_i[84][0]-1.0; 
a_i[84] [2]»2.0;  a_i [84] [3] -2 .0;  a_i [84] [4] -1 . 0 ;  a_i[8S][0]»-2.0; 
a_i[85] [2]-2.0;  a_i [85] [3] -4 .0;  a_i [85] [4] -2 . 0 ;  a_i[86][0]»-1.0; 
a_i[86]  [2]»4.0/  a_i  [86]  [4] -2 . 0;  a_i  [87]  [0] -1 . 0 ;  a_i[87][l]—l.O; 
a_i[87]  [3]— 2.0;  a_i  [88]  [0] -2 .0;  a_i  [88]  [2] -2 . 0 ;  a_i[88][3]—2.0; 
a_i[88]  [4]«1.0;  a_i  [89]  [0] -2 .0;  a_i  [89]  [2] -2 . 0 ;  a_i[89][3]=2.0; 
a_i[89]  [4]=2.0;  a_i  [90]  [0] -1 . 0 ;  a_i [90]  [3] -2 . 0 ;  a_i[90][4]-1.0; 
a_i[91]  [2]=4.0;  a_i  [91]  [3] --2 .0;  a_i [91]  [4] -2 . 0 ;  a_i[92][0]=3.0; 
a_i[92]  [2]-2.0;  a_i  [92]  [3]  — 2.0;  a_i  [92]  [4]  =2 .0;  a_i[93][0]-1.0; 
a_i[93] [2]»2.0;  a_i [93] [3] =-2 .0 ;  a_i [94] [1] =1 . 0 ;  a_i[94][2]=2.0; 
a_i[94]  [4]=1.0;  a_i[95]  [0]— 1.0;  a_i [95]  [1]  — 1 . 0 ;  a_i[95][3]-2.0; 
a_i[95]  [4]=1.0;  a_i[96]  [2]  — 2.0;  a_i  [96]  [4]  =1 . 0 ;  a_i[97][2]=2.0; 
a_i[97]  [3]=-1.0;  a_i[97]  [4]=2.0;  a_i  [98]  [1] -1 . 0 ;  a_i[98][3]=2.0; 
a_i[99]  [0]=1.0;  a_i[99]  [2]— 2.0;  a_i [99]  [3]  — 2 . 0;  a_i[100][l]—1.0; 
a_i[100]  [2]=2.0;  a_i [100]  [4] -1 . 0 ;  a_i [101]  [0]  =1 . 0 ;  a_i[101][l]=1.0; 
a_i[101]  [3]— 2.0;  a_i  [101]  [4]  =1.0;  a_i  [102]  [0]  =1 . 0 ;  a_i[102][2]—2.0; 
a_i[102] [3]=2.0;  a_i [103] [0] =2 . 0 ;  a_i [103] [3] =2 . 0 ;  a_i[104][2]=2.0; 
a_i[104]  [3]=4.0;  a_i [104]  [4] =2.0;  a_i [105]  [1] =1. 0 ;  a_i[105][3]«1.0; 
S_i[0] [0] =-171996.0;  S_i [0] [1] =-174 .2;  S_i [1] [0] =2062 . 0 ;  S_i[l][l]=.2; 
s_i[2]  [0]=46.0;  s_i[3]  [0]=11.0;  s_i [4]  [0]  =-3 . 0 ;  S_i[5][0]=-3.0; 

S_i[6]  [0]=-2.0;  S_i[7]  [0]-1.0;  S_i  [8]  [0]  =-13187 . 0;  S_i  [8]  [1]  =-1.6; 

S  i[9]  [0]=1426.0;  s_i [9]  [1] =-3 .4 ;  S_i  [10]  [0] =-517 . 0 ;  S  i [10]  [1] =1.2 ; 


s_ 

i[lll 

[0] 

=217  ; 

s_i[ll]  [l]=-.5;  a 

_i[12]  [0] =129.0 

t 

s_i[12]  [l]  =  .l; 

s_ 

i[13] 

[0] 

=48.0 

s_i[14]  [0]=-22.0 

;  S  i[15]  [0]  -17 

.0 

;  S__i[15]  [!]=-.! ; 

s_ 

i[16] 

[0] 

=  -15.0;  S_i[17]  [0]=-16. 

0;  S_i[17][l]=. 

1; 

S_i[18]  [0]=-12.0; 

s_ 

i[19] 

[0] 

=  -6.0 

s_i[20][0]=-5.0; 

S_i[21]  [0]=4.0 

• 

f 

S_i[22]  [0]=4.0; 

s_ 

i[23] 

[0] 

=  -4.0 

S_i[24][0]=1.0; 

S_i[25]  [0]=1.0; 

8 

_i[26]  [0]  =-1.0; 

s_ 

i[27] 

[0] 

=1.0; 

s_i [28]  [0] =1.0;  s 

_i[29]  [0]=-1.0; 

3 

_i[30] [0] =-2274.0 ; 

s_ 

i[30] 

[1] 

=  -  .2; 

S_i[31]  [0]  =712.0 ; 

S_i[31]  [l]  =  .l; 

S 

_i[32]  [0]=-386.0; 

s_ 

i[32] 

[1] 

=  -  .4; 

s_i[33] [0] =-301.0 

;  s_i[34]  [0]=-158 

.0;  S_i [35]  [0] =123.0 

s_ 

i[36] 

[0] 

=63.0, 

S_i[37]  [0]=63.0; 

S_i[37]  [l]=.l; 

S 

_i[38]  [0]=-58.0; 

s_ 

i[38] 

[1] 

«-.l; 

S_i[39]  [0]=-59.0; 

s_i[40]  [0]=-51 

.0 

7  S_i[41]  [0]=-38.0; 

s_ 

i[42] 

[0] 

=29.0 

S_i[43]  [0]=29.0; 

S_i[44]  [0]=-31 

.0 

;  S_i[45]  [0]=26.0; 

s_ 

.i[46] 

[0] 

=21.0 

S_i[47]  [0]=16.0; 

s_i[48]  [0]=-13 

.0 

;  S_i[49]  [0]=-10.0; 

s_ 

i[50] 

[0] 

=  -7.0 

S_i[51]  [0]=7.0; 

s  i[52]  [0]=-7.0 

! 

S_i[53]  [0]=-8.0; 

s_ 

i[54] 

[0] 

=6.0; 

s_i[55]  [0]=6.0;  s 

_i[56]  [0]  =-6.0; 

8 

_i[57]  [0]=-7.0; 

s_ 

i[58] 

[0] 

=6.0; 

S_i[59]  [0]  =-5.0; 

S_i[60]  [0]=5.0; 

S 

_i[61] [0] =-5.0; 

s_ 

i[62] 

[0] 

=-4.0, 

S_i[63]  [0]=4.0; 

S_i[64]  [0]=-4.0 

S_i[65]  [0]=-3.0; 

s_ 

i[66] 

[0] 

=3.0; 

S_i[67]  [0]  =-3.0; 

s_i[68]  [0]=-3.0 

S_i[69] [0] =-2.0; 

s_ 

i[70] 

[0] 

=-3.0, 

S_i[71]  [0]  =-3.0; 

s_i[72]  [0]  =2.0 

S_i[73]  [0]=-2.0; 

s_ 

i[74] 

[0] 

=2.0; 

S_i[75]  [0]=-2.0; 

S_i[76]  [0]=2.0; 

S 

_i[77]  [0]=2.0; 

s_ 

i[78] 

[0] 

=1.0; 

S_i[79]  [0]=-1.0; 

S_i[80]  [0]=1.0; 

S 

_i[81]  [0]=- 2.0; 

s_ 

i[82] 

[0] 

=  -1.0; 

S_i[83]  [0]=1.0; 

s_i[84]  [0]  —1.0 

} 

S  i[85]  [0]=-1.0; 

s_ 

i[86] 

[0] 

o 

II 

S_i[87]  [0]=1.0;  S 

_i[88]  [0]=1.0; 

s_ 

i[89]  [0]=-1.0; 

s_ 

i[90] 

[0] 

=-1.0; 

S_i [91]  [0] =1.0; 

S  i[92]  [0]=1.0; 

s 

_i[93]  [0]=-1.0; 

s_ 

i[94] 

[0] 

=1  -  0 ; 

S_i[95]  [0]=1.0;  S 

_i[96]  [0]=-1.0; 

s 

_i[97]  [0]=-1.0; 
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s_i[98]  tO]»-1.0;  S_i[99J  [0]«-1.0;  S_i  [100]  [0]  — 1 . 0 ;  S_itl01][0]=-1.0; 
S_i[102]  t0]=-1.0;  S_itl03]  [0] *1.0;  S_i  [104]  [0] =-l . 0 ;  S_i[105][0]=1.0; 
C_i[0]  [0]«92025.0;  C_i  [0]  [1]  =8 . 9 ;  C_i  [1]  [0]  »-895 . 0 ;  C_i[l][l]  =  .5; 

C_i[2]  [0]— 24.0;  C_i[4l  [01-1.0;  C_i  [6]  [0] -1 . 0 ;  C_i[8][0]-5736.0; 

C_i[8]  [11—3.1;  C_i[9l  [01=54.0;  C_i  [91  [1]  .  1 ;  C_i[10l[0]«224.0; 

C_i[101  [11  — .6;  C_i  [111  [0]  =-95.0;  C_i  [11]  [1]  =  .  3  ;  C_i[12][0]=-70.0; 

C_i[13]  [01=1.0;  C_i[161  [01=9.0;  C_i  [17]  [0]  =7 . 0  ;  C_i[181[0]=6.0; 

C_i[191  [01=3.0;  C_i  [20]  [01=3.0;  C_i  [21]  [0]  =-2 . 0 ;  C_i[22l[0]=-2.0; 

C_i[30]  [0]  =977.0;  C_i  [30]  [1]  =- .  5 ;  C_i  [31]  [0]  =-7 . 0 ;  C_i[32l[0]=200.0; 

C_i  [33]  [01=129.0;  c_i  [33]  [1]  =- .  1 ;  C_i  [34]  [0]  =-l .  0 ;  c_i[35][01=-53.0; 
C_i[36]  [01=-2.0;  C_i[371  [0]=-33.0;  C_i  [38]  [0]  =32 . 0 ;  C_i[391[0]=26.0; 

C_i  [40]  [01=27.0;  C_i  [41]  [0]  =16 . 0 ;  C_i  [42]  [0]  — 1 . 0 ;  C_i[431[0]=-12.0; 
C_i[44]  [0]  =13.0;  C_i[45]  [01=-1.0;  C_i  [46]  [0]  =-10 . 0 ;  C_i[47l[01=-8.0; 
C_i[481  [01=7.0;  C_i  [491  [01=5.0;  C_i  [511  [01  — 3 . 0  ;  C_i[521[01=3.0; 

C_i[53]  [01=3.0;  C_i[551  [01=-3.0;  C_i  [56]  [0]  =3 . 0  ;  C_i[57][01=3.0; 

C_i[58]  [01--3.0;  C_i  [591  [01=3.0;  C_i  [61]  [0]  =3 . 0 ;  C_i[67l[0]=1.0; 

C_i[68]  [01=1.0;  C_i[691  [01=1.0;  C_i  [70]  [0]  =1 . 0 ;  C_i[71][01=1.0; 

C_i[72]  [0] =-1.0;  C_i[73]  [01=1.0;  C_i [74]  [0] =-l . 0 ;  C_i[75l[01=1.0; 

C_i[77l  [01  =-1.0;  C_i[78]  [0]=-1.0;  C_i  [80]  [0]  =-l .  0 ;  C_i[81][01=1.0; 

C_i[83]  [01 =-1.0;  C_i[84]  [01=1.0;  C_i [851  [0] =1 . 0 ;  c_i  [88]  [0] =-l . 0 ; 
a_n[01  [41=1.0;  a__n[ll  [11=1.0;  a_n[21  [21=2.0;  a_n[21  [31  =-2.0;  a__n[21  [41=2.0 
a_n[3]  [21=2. 0;  a_n[31  [41=2.0; 

sInCOl  [01— 725.0;  s_n [0]  [1]  =224 . 0 ;  s_nll3  [0]  =523 . 0 ;  s_n [1]  [1]  =-24 . 0 ; 
s_n[2]  [01=102.0;  S_n [2]  [1]  =-47 . 0 ;  s_n[3l[01=-81.0; 

C_n  [01  [01=417.0;  C_n  [01  [11=213.0;  C_n  [1]  [0]  =61 . 0 ;  C_n[ll[ll=208.0; 
c_n[21  [01=-118.0;  c_n[2]  [1] =-41.0;  c_n [3]  [1] =32 . 0 ; 

//  Initialize  the  DE118/LE62  to  DE200/LE200  transformation  matrix************ 
m_0 [0]  [01  =  . 9999256791774783 ;  m_0 [0]  [1] =- . 0111815116768724 ; 
m_0[0] [2] =-.0048590038154553;  m_0 [1] [0] =. 0111815116959975 ; 
m_0[l]  [1]  =  .9999374845751042;  m_0  [1]  [21  —  .0000271625775175; 
m_0 [2]  [0]=. 004859003771445;  m_0 [2]  [1] =-. 000027170449221 ; 
m_0[21 [2] =.9999881946023742; 

//  Read  ECEF_WGS_84  state  vectors  from  external  files  *********************** 
if ( (j=fopen("statep.001", "r") )»=NULL) 

{ 

f close (j ) ; 

nerror ("unable  to  open  statep.OOl  in  mainO"); 

} 

if  ( (lc=fopen("statev.001",  "r")  )==NULL) 

{ 

f close (k) ; 

nerror ( "unable  to  open  statev.OOl  in  mainO"); 

} 

for ( i=0 ; i<NUMBER_OF_VECTORS ; i++ ) 

{ 

f scanf ( j , "%Lf  %Lf  %Lf  %Lf \n" , &t_0 [i] [0] , &v_0 [i] [0] , &v_0  [i]  [1] , 
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&V_0  [i]  [2]  )  ; 

fscanf  (k,  "%Lf  %Lf  %Lf  %Lf  \n“ ,  &t_0  [i]  [l],&v_0[i]  [3],&v_0[i]  [4], 

Stv__0  [i]  [5] )  ; 

t_0[i] [01 *»1000000.0;  v_0[il [01 *=1000000 . 0 ;  v_0 [il [11 *=1000000 . 0 ; 
v_0[il [21*=1000000.0; 

t_0 [il [11 *=1000.0;  v_0 [il [31 *=1000.0;  v_0 [il [4l *=1000 . 0 ; 
v_0 [il [51 *=1000.0; 

t_0[il  [01 +=522009.0;  t_0[il  [11 +=522009 . 0 ; 
if (i<5) 

{ 

printf ( "%Lf  %Lf \n" , t_0 [il  [01 , t_0 [il  [11 ) ; 

} 

if  (t_0[il  [01  !=t_0[il  [11 ) 

{ 

nerror Ctitne  element  mismatch  in  mainO"); 

} 

} 

if (fclose{j)==EOF) 

{ 

nerror ( "unable  to  close  statep.OOl  in  mainO"); 

} 

if (f close (k) ==EOP) 

{ 

nerror ( "unable  to  close  statev.OOl  in  mainO"); 

} 

//  Main  loop  to  calculate  time,  rotation  matrices,  cuid  treuisformations******** 
for  ( i=0 ; i<NUMBER_OF_VECTORS ; i++) 

{ 

bed  tdt; 

long  double  g, integer, fraction; 
long  double  t , zeta, z, theta,p [31 [31 ; 

long  doxible  l,l_prime,f  ,d,  omega, epsilon_bar,  epsilon,  c  [31  [31  ; 
long  double  utl_utc_interpolated, delta_psi=0 . 0 , delta_epsilon=0 . 0 ; 
long  double  t__u ,  h_0 ,  del ta_h ,  omega_s tar ,  lambda ,  b  [3 1  [31,  b_dot  [31  [31  ; 
long  double  xjp_interpolated,yjp_interpolated, a [31 [31 ; 
long  double  pm_0[31  [31  ,cpm_0[31  [31  ,bcpm_0[31  [31  ,b_dotcpm_0  [31  [31  ; 
long  double  abcpm_0[31 [31 , ab_dotcpm_0 [31 [31 ; 

long  double  abcpm_0_transpose  [31  [31  ,  ab__dotcpm__0_transpose  [31  [31; 

//  Convert  time  (t_0[il  [01  &  t_0[il  [ID  to  UTC/TDB/UTl  (t_0[il  [21  -  t_0[il  [51) 
//  Compute  UTC  time  in  hours  since  the  beginning  of  the  year*************** 
t_0 [i] [21 = ( (t_0 [il [01 /DAY) + (GPS_WEEK_NUMBER*7 . 0 ) + JD_GPS_WEEK_0 - 
(GPS_UTC_OFFSET/DAY) -JD_00_JAN_93 ) *24 . 0 ; 

//  Compute  Julicui  Date  in  TDT  time  rounded  to  two  decimal  places*********** 
tdt=bcd(  (  (t_0  [il  [01  /DAY)  +  (GPS_WEEK_NUMBER*7 .0)  +  JD_GPS_WEEK_0 - 

(GPS_UTC_OFFSET/DAY)  +  (TAI_UTC_OFFSET/DAY)  +  (TDT_TAI_OFFSET/DAY)  )  ,  2)  ; 
//  Compute  mean  anomaly  of  the  Earth  in  its  orbit************************** 
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g=357. 53+ ( .9856003* (real (tdt) -JD_J2000) ) ; 

//  Compute  Julian  Date  in  TDB  time***************************************** 
t_0 [i] [3] = (t_0 [i] [0] /DAY) + (GPS_WEEK_NUMBER*7 .0) + JD_GPS_WEEK_0 - 

(GPS_UTC_OFFSET/DAY) + (TAI_UTC_OFFSET/DAY) + (TDT_TAI_OFFSET/ 
DAY) +( { ( .001658*sinl (g*DEG_TO_RAD) ) + ( . 000014*sinl (2 .0*g* 
DEG_TO_RAD) ) ) /DAY) ; 

//  Compute  Julian  Date  in  UTl  time***************************************** 
if ( ( ( t_0 [i] [0] /DAY) + (GPS_WEEK_NUMBER*7 .0) +JD_GPS_WEEK_0- 
( GPS_UTC_OFFS ET / 

DAY) -2400000.5) <49248.0) 


{ 

t_0 [i] [4] = (t_0 [i] [0] /DAY) + (GPS_WEEK_NUMBER*7 .0) + JD_GPS_WEEK_0 - 
(GPS_UTC_OFFSET/DAY) + (UT1_UTC_49247/DAY) ; 

} 

if ( ( ( t_0 [i] [0] /DAY) + (GPS_WEEK_NUMBER*7 .0) + JD_GPS_WEEK_0 - 
( GPS_UTC_OFFSET/ 

DAY) -2400000.5) >=49248.0) 


// 


{ 

t_0  [i]  [4]  =  (t_0  [i]  (0]  /DAY)  +  (GPS_WEEK_NUMBER*7 . 0 )  +  JD_GPS_WEEK_0 - 
(GPS_UTC_OFFSET/DAY) + (UT1_UTC_49248/DAY) ; 

} 

Compute  Julian  Date  in  UTl  time  at  the  beginning  of  the  day************* 
fraction=modfl (t_0  ti]  [4] , &integer) ; 
if (fraction< .5) 


{ 

t_0[i] [4] =integer- .5; 

} 

if (fraction>.5) 


{ 

t_0 [i] [4] =integer+ . 5 ; 

} 

//  Compute  UTC  time  in  seconds  since  the  beginning  of  the  day************** 
fraction=modfl  (  ( (t_0  [i]  [0]  -GPS_'  ■”C_OFFSET)  /DAY)  ,  &integer)  ; 
t_0  [i]  [5]  =  (t_0  [i]  [0]  -GPS_UTC_OFir C>ET)  -  (integer*DAY)  ; 

//  Calculate  General  Precession  Rotation  Matrix  p[] []************************* 
//  Calculate  number  of  centuries  of  TDB  elapsed  since  J2000 . 0************** 
t=(t_0(i] [3] -JD_J2000) /CENTURY; 

//  Calculate  Accumulated  Precession  Angles  Adopted  by  lAU  1976************* 
zeta= ( (2306 .2181*t) + ( .30188*t*t) + ( . 017998*t*t*t ) ) *ARCSEC_TO_RAD; 
z= ( (2306 .2181*t) + (1.09468*t*t) + (. 018203* t*t*t) ) *ARCSEC_TO_RAD; 
theta=( (2004.3109*t) - ( .42665*t*t) - ( . 041833*t*t*t ) ) *ARCSEC_TO_RAD ; 
p  [0] [0]  =  (cosl (z) *cosl (theta) *cosl (zeta) ) - (sinl (z) *sinl (zeta) ) ; 
p  [0]  [1]  =  ( -cosl (z) *cosl (theta) *sinl (zeta) ) - (sinl (z) *cosl (zeta) ) ; 

P [0]  [2] =-cosl (z) *sinl (theta) ; 

p[l]  [0]  =  (sinl (z) *cosl (theta) *cosl (zeta) )  +  (cosl (z) *sinl (zeta) ) ; 
p[l] [1] = (-sinl (z) *cosl (theta) *sinl (zeta) ) + (cosl (z) *cosl (zeta) ) ; 
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p[l] [2] =-sinl (z) *sinl (theta) ;  p(2] [0] *sinl (theta) *cosl (zeta) ; 
p[2] [11 =-sinl (theta) *sinl (zeta) ;  p[2] [2] »cosl (theta) ; 

//  Calculate  Astronomical  Nutation  Rotation  Matrix  c[] []**♦***♦*************** 
//  Calculate  lAU  1980  Nutation  set  of  Fundamental  Arguments**************** 
1=  (485866. 733+ ( ( (1325 . 0*ONE_RBV) +715922 . 633 ) *t)  +  (31 . 31*t*t) +  (. 064*t*t* 
t) ) *ARCSEC_TO_RAD; 

l_prime= (1287099.804+ ( ( (99 . 0*ONE__REV) +1292581 . 244) *t) - ( .577*t*t) - ( .012* 
t*t*t) ) *ARCSEC_TO_RAD; 

f= (335778.877+ ( ( (1342 . 0*ONE_REV) +295263 . 137) *t) - (13 .257*t*t) + ( . 011*t*t* 
t) ) *ARCSEC_TO_RAD; 

d=(1072261.307+( ( (1236 . 0*ONE_REV) +1105601 . 328) *t) - 
(6.891*t*t) +( .019*t*t* 

t) ) *ARCSEC_TO_RAD; 

omega= (450160 . 28- ( ( (5 . 0*ONE_REV) +482890 . 539) *t) + (7 . 455*t*t) + ( . 008*t*t* 
t ) ) *ARCSEC_TO_RAD ; 

//  Calculate  Mean  Obliquity  of  the  Ecliptic******************************** 
epsilon_bar= (84381.448- (46.815*t) - ( . 00059*t*t) + ( . 001813*t*t*t) ) * 
ARCSEC_TO_RAD; 

//  Calculate  Angle  of  Nutation  in  Longitude******************************** 
for (n=0/n<106;n++) 

{ 

long  double  deltajpsi_i; 

delta__psi_i«  (s_i  [n]  [0]  +  (s_i[n]  [11  ♦t) )  *sinl  ( {a_i  [nj  [OJ  *1)  +  (a_i  [nj  [11* 

l_prime)  +  (a_i[n]  [2]  *f )  +  (a_i  [n]  [3]  *d)  +  (a__i  [n]  [4]*omega)); 

delta_psi+= (delta_psi_i* .0001*ARCSEC_TO_RAD) ; 

} 

//  Calculate  and  apply  Correction  to  Angle  of  Nutation  in  Longitude******** 
for (n=0;n<4;n++) 

{ 

long  double  delta_psi_c; 

delta_psi_c= (s_n[n] [0] *sinl ( (a_n[n] [0] *1) + (a_n [n] [1] *ljprime) + 

(a_n[n] [2] *f )  +  (a_n[n] [3] *d)  +  (a_n [n] [4] *omega) ) )  + 

(c_n[n]  [0]  *cosl  ( (a_n[n]  [0]  *1) +  (a_n[n]  [l]  *l_prime)  + 

(a_n[n]  [2]  *f )  +  (a_n[n]  [3]  *d)  +  (a_n[nl  [4]  *omega) ) )  ; 

delta_psi+= (delta_psi_c* . 00001*ARCSEC_TO_RAD) ; 

} 

//  Calculate  Angle  of  Nutation  in  Obliquity******************************** 
for (n=0 ;n<106 ;n++) 


long  doxible  delta_epsilon_i; 

delta_epsilon_i=  (c_i  [n]  [0]  +  (c_i  [n]  [1]  *t) )  *cosl  ( (a__i  [n]  [0]  *1)  + 
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(a_i  [n]  [1]  *l_priine)  +  (a_i  [n]  [2]  *f )  +  (a_i  [n]  [3]  *d)  + 

(a_i  [n]  [4]  *oinega) )  ; 

delta_epsilon+»  (delta_epsilon_i*  .0001*ARCSEC_TO_RAD)  ; 

} 

//  Calculate  and  apply  Correction  to  Angle  of  Nutation  in  Obliquity******* 
for (n=0;n<4fn++) 

{ 

long  double  delta_epsilon_c; 

delta_epsilon_c= (c_n[n] [1] *cosl ( {a_n[nl  [0] *1)  +  (a_n [n] [1] *l_prime) + 

{a_n(n]  [2]  *f )  +  (a_n[n]  [3]  *d) +  (a_n[n]  [4]  *oinega) ) )  + 

(s_n[n] [1] *sinl ( (a_n[n] [0] *1) +(a_n[n] [1] *l_prime) + 

(a_n[n] [2] *f ) + (a_n [n] [3] *d) + (a_n[n] [41 *omega) ) ) ; 

delta_epsilon+*  (delta_epsilon_c* .  00001*ARCSEC_TO_RAD)  ; 

} 

//  Calculate  True  Obliquity  of  the  Ecliptic******************************* 
epsilon=epsilon__bar+delta_epsilon; 

c  [0]  [0]  sscosl  (delta_psi) ;  c  [0]  [1]  »-sinl  (deltajpsi)  *cosl  (epsilon_bar) ; 
c  [0]  [2]  =-sinl  (delta__psi)  *sinl  (epsilon_bar) ; 
c [1] [0] »cosl (epsilon) *sinl (deltajpsi) ; 

c  [1]  [1]  =  (cosl (epsilon) *cosl (deltajpsi) *cosl (epsilon_bar) )  + 

(sinl  (epsilon)  *sinl  (epsilon_j3ar) )  ; 
c [1] [2] = (cosl (epsilon) *cosl (delta_psi} *sinl (epsilon_bar) ) - 
(sinl (epsilon) *cosl (epsilon_bar) ) ; 
c [2] [0] =sinl (epsilon) *sinl (delta_psi) ; 

c  [2] [1]  =  (sinl (epsilon) *cosl (delta_psi) *cosl (epsilon_bar) ) - 
(cosl  (epsilon)  *sir.l  (epsilon_bar) )  ; 
c [2] [2]  =  (sinl (epsilon) *cosl (deltajpsi)  *sinl (epsilon_bar) )  + 

(cosl (epsilon) *cosl (epsilon_bar) ) ; 

//  Calculate  Earth  Rotation  (Sidereal  Time)  Matrix  b[] []*******************♦♦ 
//  Interpolate  values  of  utl_utc****************************************** 
utl_utc_interpolated= ( ( ( ( (t_0 [i]  [2]  /24 . 0) +JD_00_JAN_93) -2400000 . 5) - 


19247.0) * (UT1_UTC_49248-UT1_UTC_49247) ) + 

//  Calculate  number  of  centuries  of  UT  elapsed 


UT1_UTC_49247 ; 
since  12h  01  JAN  2000  UTl** 


t_u=(t_0[i] [4] -2451545.0) /CENTURY; 

//  Calculate  Greenwich  Mean  Sidereal  Time  at  Oh  UTl  of  day  of  interest**** 
h_0= (24110 . 54841+ (8640184 . 812866*t_u) + ( . 093104*t_u*t_u) - ( . 0000062*t_u* 
t_u*t_u) ) *SEC_T0_RAD; 

//  Calculate  the  Equation  of  the  Ec[unoxes********************************* 
delta_h=atanl (cosl (epsilon) *tanl (deltajpsi) ) ; 

//  Calculate  Earth  rotation  rate  in  a  precessing  reference  frame********** 
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oinega_star=OMEGA_PRIMB+.000000000007086+(  .0000000000000043*t_u)  ; 

//  Calculate  Longitude  of  the  Zero  Meridiem  from  the  true  vernal  equinox**^ 
lambda=h_0+delta_h+ (omega_star* (t_0 [i]  [5] +utl_utc_interpclated) ) ; 
b[0]  [0]  scosl  (lambda)  ;  b[0]  [1]  ssinl  (lambda)  ;  b[0]  [21=0.0; 
b[l]  [0]  s-sinl  (lambda)  ;  b  [1]  [1]  =cosl  (lambda)  ;  b[l][2]=0.0;  b[2][0]=0.0; 
b[21  [11=0.0;  b[21  [21=1.0; 

//  Calculate  Change  in  Earth  Rotation  (Sidereal  Time)  Matrix  b_dot[l  [l********| 
b_dot[01  [01 =-omega_star*sinl (lambda) ; 
b_dot [01 [11  =omega_star*cosl (lambda) ;  b_dot [01 [21=0.0; 
b_dot [11 [01 =-omega_star*cosl (lambda) ; 

b_dot [11 [11 =-omega_star*sinl (lambda) ;  b_dot [11 [21 =0 . 0 ;  b_dot [21 [01 =0 . 0 ; 
b_dot[21  [11=0.0;  b_dot[2l  [21=0.0; 

//  Calculate  Polar  Motion  Rotation  Matrix  a[l  [l  *******************************j 
//  Interpolate  values  of  x__p  &  y_p  emd  convert  to  radiams****************** 
x_p_interpolated= ( ( ( ( ( (t_0[il  [2] /24 . 0) +JD_00_JAN_93) -2400000.5) - 


49247.0) * (X_P_49248-X_P_49247) ) +X_P_49247) * 

ARCSEC_TO_RAD; 

y_p_interpolated= ( ( ( ( ( (t_0[i) [21 /24 . 0) +JD_00_JAN_93) -2400000.5) - 


49247.0) ♦(¥  P  49248-Y  P  49247) )+Y  P  49247)* 


ARCSEC  TO  RAD; 


a  [01  [01  =cosl  (x_j)__interpolated)  ; 

a [01 [11  =sinl (x_p_interpolated) *sinl (y_p_interpolated) ; 
a  [01  [21  =sinl  (x_p__interpolated)  *cosl  (y_p_interpolated)  ;  a  [11  [01  =0 . 0  ; 
a [11 [11 =cosl (y_p_interpolated) ;  a [11  [2]  =-sinl (y_p_interpolated) ; 
a  [21  [01 =-sinl (xjp_interpolated) ; 

a [21 [11 =cosl (x_p_interpolated) ♦sinl (yjp_interpolated) ; 
a [21 [21 =cosl (x_p_interpolated) *cosl (y_p_interpolated) ; 

//  Perform  matrix  multiplication  to  form  transformation  matrices************** 
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pm_0 [2]  [01 )  / 

cpm_0  [1]  [1]  »(c[l]  [01  *pm_0  [01  [11 )  +  (c[ll  [11  *pm_0  [11  [11 )  +  (c[ll  [21  * 
pm_0  [21  [11 )  ; 

cpm_0[ll  [21  =  (c[ll  [01*pm_0[01  [21)  +  (c[ll  [ll*pm_0[ll  [21)  +  (c[ll  [21* 
pm_0[21  [21 )  ; 

cpm_0  [21  [01  =  (c  [21  [01  ♦pm_0  [01  [01 )  +  (c  [21  [11  *pm_0  [11  [01 )  +  (c  [21  [21  * 
pm_0[21  [01)  ; 

cpm_0[21  [ll  =  (c[21  [01*pin_0[01  [ll)  +  (c[21  [ll*pm_0[ll  [11  )  +  (c[21  [21* 
pm_0  [21  [11 )  ; 

cpm_0[21  [21  =  {c[21  [01*pm_0[01  [21)  +  (c[21  [ll*pm_0[ll  [21)  +  (c[21  [21* 
pm_0 [21  [21 )  ; 

bcpin_0[01  [01  =  (b[01  [01*cpm_0[01  [01)  +  (b[0l  [11  *cpm_0  [11  [01  )  +  (b[01  [21* 
cpm_0  [21  [01 )  ; 

bcpm_0[01  [ll  =  (b[01  [01*cpm_0[01  [ll)  +  (b[0l  [11  *cpm_0  [11  [ll)  +  (b[01  [21* 
cpm_0 [21  [11 )  ; 

bcpin_0[01  [21  =  (b[01  [01*cpm_0[01  [21)  +  (b[01  [ll*cpm_0[ll  [21)  +  (b[01  [21* 
cpin_0  [21  [21 )  ; 

bcpm_0[ll  [01  =  (b[ll  [01*cpm_0[01  [01)  +  (b[ll  [ll*cpm_0[ll  [01)  +  (b[ll  [21* 
cpm_0 [21  [01 )  ; 

bcpm_0[ll  [ll  =  (b[ll  [01*cpm_0[01  [ll)  +  (b[ll  [11  *cpra_0  [11  [ll)  +  (b[ll  [21* 
cpm_0  [21  [11 )  ; 

bcpm_0[ll  [21  =  (b[ll  [01*cpm_0[01  [21)  +  (b[ll  [11  *cpm_0  [11  [21)  +  {b[ll  [21* 
cpm_0  [21  [21 )  ; 

bcpm_0[21  [01  =  (b[21  [01  *cpm_0  [01  [01)  +  (b[21  [ll*cpm_0[ll  [01)  +  (b[21  [21* 
cpin_0  [21  [01 )  ; 

bcpm_0[21  [ll  =  (b[2l  [01  *cpm_0  [01  [ll)  +  (b[2l  [ll*cpm_0[ll  [11  )  +  (b[21  [21* 
cpm_0 [21  [11 )  ; 

bcpm_0[21  [21  =  (b[21  [01  *cpm_0  [01  [21  )  +  (b[21  [11  *cpTn_0  [11  [21  )  +  (b[21  [21* 
cpm_0 [21  [21 )  ; 

b_dotcpm_0 [01  [01  =  (b_dot [01  [01 *cpm_0  [01  [01 )  +  (b_dot  [01  [11  *cpm_0 [11  [01 )  + 
(b_dot [01 [21 *cpm_0 [21 [01 ) ; 

b_dotcpin_0  [01  [11  =  {b_dot  [01  [01  *cpm_0  [01  [11  )  +  {b_dot  [01  [11  *cpm_0  [11  [11 )  + 
(b_dot  [01  [21  *cpm_0  [21  [11 )  ; 

b_dotcptn_0  [01  [21  =  (b_dot  [01  [01  *cpm_0  [01  [21 )  +  (b_dot  [01  [11  *cpm_0  [11  [21 )  + 
(b_dot [01  [21 *cpm_0 [21  [21 )  ; 

b_dotcpTn_0  [11  [01  =  (b_dot  [11  [01  *cpm_0  [01  [01 )  +  (b_dot  [11  [11  *cpm_0  [11  [01 )  + 
(b_c?ot  [11  [21  *cpm_0  [21  [01 )  ; 

b_dotcptn_0  [11  [11  =  (b_dot  [11  [01  *cpm_0  [01  [11 )  +  (b_dot  [11  [11  *cpm_0  [11  [11  )  + 
(b_dot  [11  [21  *cpm_0  [21  [11 )  ; 

b_dotcpm_0  [11  [21  =  (b_dot  [11  [01  *cptn_0  [01  [21 )  +  (b_dot  [11  [11  *cpm_0  [11  [21  )  + 
(b_dot [11  [21  *cpm_0 [21  [21 )  ; 

b_dotcpm_0 [21  [01  =  (b_dot [21  [01 *cpm_0 [01  [01 )  +  (b_dot [21  [11  *cpm_0 [11  [01 )  + 
(b_dot [21  [21 *cpm_0 [21  [01 )  ; 

b_dotcpm_0 [21  [11  =  (b_dot [21  [0] *cpm_0  [01  [11 )  +  (b_dot [21  [11  *cpm_0 [11  [11 )  + 
(b_dot [21  [21  *cpm_0 [21 [ll ) ; 

b_dotcpm_0 [21  [21  =  (b_dot [21  [01 *cpm_0 [01  [21 )  +  (b_dot [21  [ll *cpm_0 [11  [21 )  + 


77 


(b_dot  [21  [21  *cpm_0  [21  [21 )  ; 

cd3cpni_0[01  [01  =  (a  [01  [01  *bcpm_0  [01  [01)  +  (a[01  [11  *bcpm_0  [11  [01)  +  (a[01  [21* 
bcpm_0 [21 [01 ) ; 

abcpm_0  [01  [11  =  (a  [01  [01  *bcpm_0  [01  [11 )  +  (a  [01  [11  *bcpm_0  [11  [11 )  +  (a[01  [21  * 
bcpm_0 [21 [11 ) ; 

abcpm_0(01  [21  =  (a[0]  [01  *bcpm_0  [01  [21)  +  (a[01  [11  *bcpm_0  [11  [21)  +  (a[01  [21* 
bcpm_0 [21  [21 ) ; 

abcpm__0[ll  [01  =  (a [11  [01  *bcpin_0  [01  [01)  +  (a[ll  [11  *bcpm_0  [11  [01)  +  (a[ll  [21* 
bcpm_0  [21  [01 )  ; 

e03cpm_0[ll  [11  =  (a  [11  [01  *bcpm_0  [01  [ll)  +  (a[ll  [1]  *bcpm_0  [11  [ll)  +  (a[ll  [21* 
bcpm_0 [21 [11 )  ; 

abcpm_0  [11  [21  =  (a  [11  [01  *bcpm_0  [01  [21 )  +  (a  [11  [11  *bcpm_0  [11  [21 )  +  (a[ll  [21  * 
bcpm_0 [21 [21 ) ; 

cU3cpra_0  [21  [01  =  (a  [21  [01  *bcpm_0  [01  [01 )  +  (a  [21  [11  *bcpm_0  [11  [01 )  +  (a[21  [21  * 
bcpm_0  [21  [01 )  ; 

c0Dcpm__0  [21  [11  =  (a  [21  [01  *bcpm_0  [01  [ll)  +  (a[21  [11  *bcpm_0  [11  [ll)  +  (a[21  [21* 
bcpm_0 [21  [11 ) ; 

abcpm_0[21  [21  =  (a  [21  [01  *bcpm_0  [01  [21)  +  (a[21  [11  *bcpm_0  [11  [21  )  +  (a[21  [21* 
bcptn_0  [21  [21 )  ; 

cQ3_dotcpm_0  [01  [01  =  (a  [01  [01  *b_dotcpm_0  [01  [01)  +  (a[0l  [11  *b__dotcpm_0  [11  [01  )  + 

(a[01  [21  *b_dotcpm_0[21  [01 )  ; 

eUD__dotcpm_0  [01  [11  =  (a  [01  [0]  *b_dotcpm_0  [01  [l])-t-(a[01  [11  *b_dotcpm_0  [11  [11 )  + 

(a[0] [21 *b_dotcpra_0 [21 [ll ) ; 

ab__dotcpin__0  [01  [21  =  (a  [01  [01  *b__dotcpin_0  [01  [21)"»-(a[01  [11  *b_dotcpra_0  [11  [21 )  + 

(a[01  [21  *b__dotcpm_0 [21  [21 )  ; 

ed)_dotcpm_0  [11  [01  =  (a  [11  [01  *b_dotcpra_0  [01  [0] )  +  (a  [ll  [11  *b_dotcpm_0  [11  [01 )  + 

(a  [11  [21  *b__dotcpm_0  [21  [01 )  / 

ab_dotcpm_0 [ll  [ll  = (a [ll  [01 *b_dotcpm_0 [01  [1] )  +  (a [11  [11  *b_dotcpm_0 [ll [ll )  + 

(a  [11  [21  *b__dotcpm_0  [21  [11 )  ; 

ab_dotcpm_0  [ll  [21  =  (a  [11  [01  *b_dotcp!n_0  [01  [21)  +  (a[ll  [ll  *b_dotcpm_0  [11  [21  )  + 

(a [11 [21 *b_dotcpm_0 [21 [21 ) ; 

ab_dotcpm_0 [21  [01  =  (a [21  [01 *b_dotcpm_0 [01 [01 )  +  (a [21  [ll *b_dotcpm_0 [ll  [01 )  + 

(a [21  [21 *b_dotcpm_0  [21  [01 ) ; 

ab_dotcpm_0  [21  [11  =  (a  [21  [01  *b_dotcpm_0  [01  [l])  +  (a[2l  [11  *b_dotcpm_0  [ll  [ll)  + 

(a  [21  [21  *b_dotcpra_0  [21  [11 )  ; 

ab_dotcpm_0  [21  [21  =  (a  [21  [01  *b_dotcpra_0  [01  [21 )  +  (a  [21  [ll  *b_dotcpm_0  [11  [21 )  + 

(a [21  [21  *b_dotcpm_0[21  [21 ) ; 

//  Transpose  the  transformation  niatrices************************************** 
abcpm_0_transpose [01 [01 =abcpm_0 [01 [0] ; 
abcpm_0_transpose[01 [1] =abcpm_0 [ll [01 ; 
cdjcpm_0_treuispose  [01  [2]  =abcpm_0  [21  [01  ; 
abcpm_0_transpose  [11  [0]  =cd3cpm_0  [01  [ll  ; 
abcpm_0_transpose[ll [11 =abcpra_0 [ll [11 ; 
abcpm_0_trcuispose  [ll  [21  =abcpm_0  [21  [11  ; 
abcpm_0_treuispose  [21  [01  =abcpm_0  [01  [21  ; 
abcpm_0_transpose  [21  [ll  =eUDCpra_0  [11  [21  ; 
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2Q3cpm_0_tran8pose[2]  [2]  >abcpai_0  [2]  [2]; 
ab_dotcpm_0_tran8pose to] [0] -ab^dotcpm^O [0] [O] ; 
ab_dotcpm_0_tr2ui8po8e  [0]  (11  «ab_dotcpm_0  [l]  [0]  ; 
edb_dotcpm_0_tran8po8e[0]  (2]  ■ab__dotqpm_0  [2]  [0]; 
ab_dotcpra_0_tran8po8e  [1]  (01  ■ab_dotcjan_0  (01  (11  ; 
ab_dotcpin_0_tran8po8e  (1]  (11  »ab__dotcpm__0  (11  (1]  ; 
cd)_dotcpm_0_tran8po8e  (1]  (21  «ab_dotcpm_0  (21  (1]  ; 
ab_dotcpm_0_traui8po8e(2]  (0]  ■ab__dotcpm_0  (0]  (21  ; 
ab_dotcpm__0_tran8po8e(21  (11  »ab_dotcpm_0{ll  (21  ; 
ab_dotcpm_0_tran8pose(2]  (21  »ad)_dotcpin_0  (2]  (21  ; 

//  Perform  BCEF_WGS_84  to  ECI_M50  po8ition  coordinate  traui8formation*********i^ 
v_l  (il  (01  =  (abcpm_0_tran8po8e  (01  (01  *v_0  (il  (01 )  +  (aQ3cpm_0_tram8po8e  (0]  (11  * 
v_0  (il  (11 )  +  (abcpm_0_tran8po8e(01  (21  *v__0  (il  (21 )  ; 
v_l  (il  (1]  =  (ad)cpm_0_tran8po8e  (1]  (01  *v_0  (il  (01 )  +  (abcpm_0_tram8po8e  (11  (1]  * 
v_0 (il (11 ) + (abcpm_0_tran8po8e (11 (21 *y_0 (il (21 ) ; 
v_l  (il  (21  =  (abcpm_0_tran8po8e  (2]  (01  *v_o  (il  (01 )  +  (cd)Cpm_0_tran8po8e  (21  (l]  ♦ 
v_0  (il  (11 )  +  (abcpm_0_tran8pO8e(2]  (21  *v_0  (il  (21 )  ; 

//  Perform  ECEF_WGS_84  to  ECI_M50  velocity  coordinate  traneformation* ********* 
v_l (il (31 = (ab_dotcpm_0_tran8po8e (01 (01 *v_0 (il (01 ) + 
(ab_dotcpm_0__traui8po8e  (01  (11  *v__0  (il  (11 )  + 
(ab_dotcpm__0_tranapo8e  (01  (21  *v__0  (il  (21 )  + 

(abcpm_0_tr8ui8po8e  (01  (01  *v_0  (il  (31 )  +  (abcpm_0__tran8po8e  (0]  (1]  * 
v_0  (il (4] )  +  (abcpm_0_tran8po8e (0] (21  *v_0 (il (51 ) ; 
v__l  (il  (41  =  (ab_dotcpm_0__tran8pose  (ll  (01  *v_0  (il  (01 )  + 
(aU3_dotcpm_0_tran8po8e  (11  (11  *v_0  (il  (ll )  * 
(cUD_dotcpm_0_tran8po8e  (11  (21  *v_0  (i]  (21 )  + 

(abcpm_0_tram8po8e  (11  (01  *v_0  (i]  (3] )  +  (ed>cpm_0_treuispo8e  (1]  (1]  * 
v_0  (il  (41 )  +  (cd)cpm_0_tran8po8e  (11  (2]  *v_0  (il  (5] )  ; 
v_l (i]  (51  = (ab_dotcpm_0_tran8po8e (21  (01  *v_0 (i]  (01)  + 
(ed3_dotcpm_0_tremspose  (2]  (11  *v_0  (il  (1] )  + 
(ab_dotcpm_0_transpose (2] (2] *v_0 (il (21 ) + 

(ed)cpm_0_transpose  (21  (01  *v_0  (il  (31 )  +  (abcpm_0_transpose  (2]  (1]  * 
v_0  (il  (41 )  +  (abcpm_0_trauispose  (21  (2]  *v_0  (il  (51 )  ; 

//  Write  status  information  to  an  external  file 
if ( { j=fopen("statout . 001" , "a" ) ) ==NULL) 

{ 

f close (j ) ; 

nerror ( "xinadale  to  open  statout.OOl  in  mainO"); 

} 

if (i<l) 

{ 

long  double  teotl , te8t2 , test! , test4 , test5 , tests , test? , test8 , test9 ; 

long  doxible  testlO ,  testll ,  testl2 ,  testl3 ,  testl4 ,  testlS ,  testlS ,  testl?  ; 

testl=xjp_interpolated/ARCSEC_TO_RAD,* 

test 2  =y_p_interpolat ed/ARCSEC_TO_RAD ; 

tes t3 = - lambda/DEG_TO_RAD ; 
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t  es  1 4  »  -  eps  ilon__bar /DEG_TO__RAD  ; 
testS-delta_psi/ARCSBC__TO_RAD; 
teste ■epsilon/DBG_TO_RAD ; 
test7»zeta/ARCSEC_T0_RAD ; 
test 8  a - theta/ARCSEC_TO_RAD ; 
test9»z/ARCSEC_T0__RAD  ; 

fprintf (j , "Rotate  about  the  Y1  eucis  by  %  16.13Lf\".\n",testl); 

fprintf (j , "Rotate  about  the  X2  eucis  by  %  16 . 13Lf\" . \n" , test2) ; 

fprintf (j , "Rotate  about  the  Z3  eucis  by  %  16.13Lf  degrees.\n",test3); 

fprintf  (j ,  "Rotate  ed)out  the  X4  eucis  by  %  16.13Lf  degrees.\n",test4); 

fprintf (j , "Rotate  about  the  Z5  axis  by  %  16. 13Lf \" .\n", tests) ; 
fprintf (j , "Rotate  about  the  X6  eucis  by  %  16.13Lf  degrees.\n",test6); 
fprintf (j , "Rotate  about  the  Z7  eucis  by  %  16.13Lf\".\n",test7); 

fprintf (j , "Rotate  about  the  Y8  axis  by  %  16.13Lf\" .\n", test8} ; 

fprintf (j , "Rotate  about  the  Z9  eucis  by  %  16.13Lf\".\n",test9); 

fprintf (j , "Rotate  about  the  ZIO  axis  by  -0.320288870000  degrees.\n"); 

fprintf  (j ,  "Rotate  ed)out  the  Yll  eucis  by  0.278405860000  degrees. \n" )  ; 

fprintf  (j ,  "Rotate  ediout  the  Z12  eucis  by  -0.320381690000  degrees . \n" )  ; 

} 

if (f close ( j ) ==EOF) 

{ 

nerror ("unable  to  close  statout.OOl  in  roainO"); 

} 

} 

//  Write  ECI_M50  state  vectors  to  an  external  file**************************** 
if  ( ( j *f open ("stateout. 001",  "w")  }«>:NULL) 

{ 

fclose( j) ; 

nerror ( "unable  to  open  stateout.OOl  in  mainO"); 

} 

for  ( i=0 ; i<NOMBER_OF_VECTORS ; i++ ) 

{ 

fprintf (j, "%Lf  %  .12Le  %  .12Le  %  .12Le\n  %  .12Le  %  .12Le  %.12Le\n", 
t_0[i]  [2]  ,v_l[i]  [0]  ,v_l[i]  [1]  ,v_l[i]  [2]  ,v_l[il  [3]  ,v_l(i]  [4]  ,v_l[i]  [51 )  ; 

} 

if (fclose( j) =*EOF) 

{ 

nerror ( "unable  to  close  stateout.OOl  in  mainO"); 

} 

//  Write  ECEF_WGS_84  state  vectors  in  Kalman  Filter  input  file  format********* 
if ( (j=fopen("sv0.001", "w") ) =*NULL) 

{ 

fclose(j) ; 

nerror ( "uned>le  to  open  svO.OOl  in  mainO"); 

} 
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for  ( i»0 ;  i<NUMBBR_OP_VBCTORS ;  i++ ) 

{ 

int  minutes; 

long  double  seconds; 

minutes> (int) (t_0(i] [51/60.0); 

seconds=t_0 [i] [5] - (60 . 0* (float)minutes) ; 

minutes-=( (int) (t_0[i] [5] /3600.0)*60) ; 

if(i>0) 

{ 

if ( (t_0 [i] [0] -t_0 [i-1] [0] ) >longest) 

{ 

longest=t_0 [i] [0] -t_0 [i-1] [0] ; 
print f ( "\n%Lf " , longest) ; 

} 

} 

fprintf (j, "%d  %S.2Lf  %  15.6Lf  %  15.6Lf  %  15.6Lf  %  12.6Lf  %  12.6Lf  % 
12.6Lf\n''. 

minutes,  seconds,  v_0  [i]  [0l,v_0[i]  [l],v_0[i]  [2],v_0[i]  [31, 
v_0[i]  [41  ,v_0[il  [51)  ; 

} 

if (fclose(j)=»EOP) 

{ 

nerror ("unable  to  close  svO.OOl  in  mainO"); 

} 

} 
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/*****************«*************************************************«*«** ! 

/*  */ 

I*  This  program  compares  state  vectors  and  computes  the  magnitude  of  */ 

/*  their  position  and  velocity  difference.  */ 

/♦  */ 

{*  Author:  LT  Stephen  P.  Rehwald,  USN  */ 

/*  Date:  20  March  1994  */ 

/*  */ 

/*  Functions  from  "Numerical  Recipes  in  C",  Press,  W.H.,  et  al,  */ 

/*  Cambridge  University  Press,  1988  were  utilized  in  this  program.  */ 

/*  */ 

/********★************************************************«************** ^ 

/************************************************************************ ! 

/*  NUMBER__OF_VBCTORS  is  the  nimiber  of  position/velocity  vectors  to  ♦/ 

/*  process  from  the  files  used  as  input  to  this  program.  Each  file  */ 

/*  must,  at  a  minimum,  contain  this  number  of  vectors.  STATBOUT.OOi  an^^  */ 
/*  REFTRAJ.OOl  specify  the  names  of  the  input  files.  M50DIPFS.001  */ 

/*  specifies  the  name  of  the  output  file.  ♦/ 


/****«*****************************♦*■************************************ y 

#define  NUMBER_OF_VECTORS  139 

#include  <math.h> 

#include  <stdio.h> 

#include  <stdlib.h> 

void  nerror (char*) ; 

FILE  *j,*lc; 
int  i; 

long  double  t_0  [NUMBER_OF_VECTORS]  [2]  ,  v_dif  f  [NUMBER_OF_VECTORS]  [2]  ; 
long  double  v_0  (NUMBER_OF_VECTORS]  [6J ,  v_l  [NDMBER_OP_VECTORS]  [6]  ; 

/*  */ 

/*  nerror  is  Numerical  Recipes  standard  error  handler*/ 

/*  */ 


void  nerror (char  error_text [] ) 

{ 

fprintf (stderr, "Numerical  Recipes  run-time  error. . .\n")  ; 
f print f  (stderr,  "%s\n" ,  error_text)  ,- 
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fprintf (stderr, " . . .now  exiting  to  system _ \n")  ; 

exit (1) ; 


void  main (void) 

{ 

//  Read  GPS_M50  cuad  Reference  Trajectory  state  vectors  from  external  files**** 
if ( {j=fopen("stateout.001", "r“) )=«NULL) 

{ 

f close (j ) ; 

nerror ( "unable  to  open  stateout.OOl  in  mainO"); 

} 

if ( (k=fopen("reftraj .001", "r") )=»NULL) 

{ 

f close (k) ; 

nerror ("unodjle  to  open  reftraj.OOl  in  mainO"); 

} 

for ( i=0 ; i<NUMBER_OF_VECTORS ; i++ ) 

{ 

fscanf ( j , "%Lf  %Le  %Le  %Le  %Le  %Le  %Le\n" , &t_0 (i]  t0],&v_0[i]  [0], 

&V_0  [i]  [11  ,  &v_0  [i]  [2]  ,  &v_0  [i]  [31  ,  &v_0  [il  [41  ,  &V_0  [il  [51  )  ; 
fscanf  (k,  "%Lf  %Le  %Le  %Le  %Le  %Le  %Le\n",  &t__0  [il  [ll,&v_l[il  [01, 

&v_l  [il  [11  ,  &v_l  [il  [21  ,  &v_l  [il  [31 ,  &v_l  [il  [41 ,  &v_l  [il  [51 )  ; 

//if (i>0) 

{ 

if  (t_0[il  [01  !=t_0[il  [ID 

{ 

nerror  ("time  element  mismatch  in  mainO"); 

} 

} 

} 

if (fclose ( j ) =*EOF) 

{ 

nerror ( "unable  to  close  stateout.OOl  in  mainO"),* 

} 

if ( fclose (k) ==EOF) 

{ 

nerror  ("unable  to  close  reftraj.OOl  in  mainO"),* 

} 

//  Main  loop  to  calculate  time,  rotation  matrices,  and  transformations******** 
for  ( i=0  ,*  i<NUMBER_OF_VECTORS  ,*  i+*i* ) 

{ 

//  Compute  the  magnitude  of  the  position  difference*************************** 
v_dif f [il  [01 =sqrtl (powl ( (v_l [i]  [01 -v_0  [il  [01  ) , 2 . 0) +powl ( (v_l  [il  [11  - 

v_0  [i]  [11  )  ,  2 .0)  +powl  (  (v_l  [il  [21  - 

V  0[il  [21 )  ,2.0) )  ; 
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//  Confute  the  magnitude  of  the  velocity  difference*************************** 
v_diff  [i]  [1]  ■sqrtl  (powl  ( (v_ltil  [3]  -v_0[il  [3] )  ,2.0)+powl{  (v_l[i]  [4]  - 

v_0  [il  t4j )  ,2.0)  +powl  ( [i]  [5]  - 

v_0[il  [5]),2.0))  ; 

} 

//  Write  differences  to  an  external  file************************************** 
if { (j=fopen("m50diffs.001". "w") )»«NULL) 

{ 

fclose(j) ; 

nerror  ( "unable  to  open  DIFFERENCB__FILE  in  main  ( )  " )  ; 

} 

for ( i=0 ; i<NUMBER_OF_VBCTORS ; i++) 

{ 

f print f (j, "%Lf  %l6.13Lf 

%16.13Lf\n",t_0[i]  [01  , v_diff  [i]  [0]  , v__dif f  [i]  [1] )  ; 

} 

if (fclose(j)==EOF) 

{ 

nerror  ( "iinable  to  close  DIFFEREMCE_FILE  in  main  ( )  " )  ; 

} 
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APPENDIX  C 


The  plots  appearing  in  this  appendix  depict  the  number  of  usable  GPS  satellites 
being  tracked  by  the  recevier.  If  more  than  three  usable  satellites  are  being  tracked,  the 
indication  is  given  either  that  the  receiver  is  doing  position  fixes,  or  that  PDOP  is  too 
high.  VDOP,  HDOP,  and  POOP  are  also  shown  (PDOP  is  a  combination  of  VDOP  and 
HOOP)  at  the  corresponding  times.  The  segment  at  the  center  (right  to  left)  of  each  plot 
corresponds  to  the  data  chosen  for  analysis  in  Chapter  III,  Table  3.  For  comparison  pur¬ 
poses  the  times  shown  in  Figures  S-10  equate  to  times  ai^iearing  on  the  TANSGRAPH 
plots  as  follows: 

•  6261.434514  hours  is  FRI:2 1:26:04.25 

•  6261.507708  hours  is  FRI:2 1:30:27.75 

•  6261.893542  hours  is  FRI:21:53:36.75 

•  6261.988264  hours  is  FRI:2 1:59: 17.75 

•  6265.089236  hours  is  SAT:01:05:21.25 

•  6265.193264  hours  is  SAT:01:11:35.75 
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Figure  C-4:  File  P.005  TANSGRAPH  Plot  Showing  PDOP 
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Figure  C-6:  File  P.008  TANSGRAPH  Plot  Showing  PDOP 
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APPENDIX  D 


This  Kalman  Filter  is  designed  for  a  user-specified  number  of  data  points  which,  in 
this  case,  is  found  in  a  file  named,  fltiest2.dat  in  ASCII  format.  Flttest2.dat  is  included 
in  this  appendix.  It  can  be  found  beginnii^  on  the  third  page  of  Af^ndix  D.  The  func¬ 
tion,  getvalsO.  reads  the  flttest2.dat  into  the  af^ropriate  vectors:  tsec,  X,  XD,  Y,  YD, 
Z,  and  ZD.  This  short  function  is  included  at  the  end  of  the  filter  program. 


load  flttest2.dat 

[tsec,  X,  XD,  Y,  YD,  Z,  ZD]  =  gctvate(  i,  flttest2); 

^Initialize  variables 

A=[01;00]; 

B=[0  ir; 

I=cye(2); 

C«[l  0;0  1]; 

Tf=4; 

dt=l; 

[Phi,Del]=c2d(A,B,dt); 

PIckml  =  le6*I; 

R=I; 

Q=.01*I; 
kmax=72  ; 
u=zeros(2,kmax); 

xkk=zeros(2,kniax+ 1); 
xkkml  =zeros(2,kmax+ 1); 
ykk=zeros(2,kmax+ 1); 
ykkml  =zeros(2,kmax+ 1); 
zkk=zeros(2,kmax+ 1); 
zkkml  =zeros(2,kmax+ 1); 

xhat=:zeros(2,kmax-l- 1); 
yhat=zeros(2,kmax+ 1); 
zhat=zeros(2,kmax+ 1); 
time  -  zeros(  1 ,  kmax) ; 
t=zeros(l,kmax); 
g=zeros(4,kmax); 
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%Bi^in  the  fUtering  process 
for  (i»l:kinax) 

Phi(l,2)*ts^i+ l)-tsec(i); 

G=Pkkml*C’*i!iv(C*Pldanl*C’  +R); 
Pkk=a-G*C)*Pkkml; 

Pkkinl=Phi*Pkk*Phi’  +  Q; 

xkkO.i)  =  xkkml(:,i)  +G*([X(i,l)  XD(i,l)r-C*xkkml(:,i)); 
xkkinl(:,i+ l)=Phi*xkk(:.i); 

ylck(:.i)=  ykkinl(:.i)  +G*([Y(i,l)  YD(i.l)]’-C*ykkml(:,i)); 
ykkinl(:,i+ l)=PhiSkk:(:,0; 

zkk(;,i)=  zkkml(:,i)  +G*([Z(i,l)  ZD(i.l)l’-C*zkkml(:,i)); 
zkkml(:  ,i+ 1) =Phi*zkk(:,i); 

xhat(:,i)  =  C*xkk(:,i); 
yhat(:,i)  =  C*ykk(:,i); 
zhat(:,i)  =  C*^(:,0; 

t(i)=time(i); 

tiiiie(i+l)*  time(i)  +  dt; 
g(l,i)=G(l,l); 
g(2,i)=Ga.i); 
g(3.i)=G(1.2); 
g(4.i)=G(2,2); 
end; 


%This  getvalsQ  function  woilcs  for  98  elements 

function  [tsec,  X,  XD,  Y,  YD,  Z,  ZD]  =  getvalsf  j ,  fittest) 

i  =1:98; 

tsec  =  60*  flttest2(i,  1)  +  flttest2(l,  2); 

X  =  flttest2(i,  3); 

XD  =  flttcst2(i,  6); 

Y  =  flttest2(i,  4); 

YD  =  flttest2(i.  7); 

Z  =  flttest2(i,  5); 

ZD  =  flttes^Ci,  8); 


93 


X 

X 

26 

04.250 

-6671.363 

276.288 

196.336 

-0.384 

-6.293 

-3.675 

26 

06.750 

-6672.305 

260.552 

187.146 

-0.364 

-6.294 

-3.676 

26 

09.250 

-6673.186 

244.821 

177.958 

-0.344 

-6.294 

-3.676 

26 

11.750 

-6674.019 

229.087 

168.770 

-0.324 

-6.295 

-3.677 

26 

14.250 

-6674.797 

213.356 

159.582 

-0.304 

-6.296 

-3.677 

26 

16.750 

-6675.528 

197.622 

150.391 

-0.284 

-6.296 

-3.678 

26 

19.750 

-6676.331 

178.744 

139.364 

-0.240 

-6.297 

-3.678 

26 

22.250 

-6676.958 

163.002 

130.168 

-0.240 

-6.297 

-3.679 

26 

24.750 

-6677.533 

147.260 

120.970 

-0.220 

-6.298 

-3.679 

26 

29.250 

-6678.442 

118.921 

104.413 

-0.184 

-6.298 

-3.680 

26 

33.750 

-6679.185 

90.579 

87.855 

-0.148 

-6.299 

-3.680 

26 

36.250 

-6679.531 

74.833 

78.653 

-0.128 

-6.299 

-3.680 

26 

38.750 

-6679.821 

59.087 

69.453 

-0.108 

-6.299 

-3.f31 

26 

40.750 

-6680.022 

46.489 

62.092 

-0.092 

-6.299 

•3.681 

26 

43.250 

-6680.238 

30.737 

52.887 

-0.072 

-6.299 

-3.681 

26 

45.250 

•6680.367 

18.138 

45.524 

-0.055 

-6.299 

•3.681 

26 

47.750 

•6680.485 

2.389 

36.320 

-0.035 

-6.299 

-3.681 

26 

49.750 

-6680.535 

-10.207 

28.960 

-0.019 

•6.299 

-3.681 

26 

52.250 

-6680.556 

-25.952 

19.758 

0.001 

-6.299 

-3.681 

26 

54.250 

-6680.532 

-38.549 

12.398 

0.017 

-6.299 

-3.681 

26 

57.750 

-6680.425 

-60.594 

-0.486 

0.045 

-6.298 

-3.681 

27 

00.250 

-6680.295 

-76.343 

-9.692 

0.065 

-6.298 

-3.681 

27 

06.250 

-6679.761 

-114.129 

-31.778 

0.113 

-6.297 

-3.681 

27 

09.250 

-6679.394 

-133.025 

-42.823 

0.137 

-6.297 

-3.681 

27 

11.750 

-6679.024 

-148.767 

-52.023 

0.157 

-6.296 

-3.680 

27 

13.750 

-6678.699 

-161.361 

-59.384 

0.173 

-6.296 

-3.680 

27 

16.250 

-6678.238 

-177.096 

-68.585 

0.193 

-6.295 

-3.680 

27 

18.750 

-6677.734 

-192.836 

-77.787 

0.213 

-6.295 

-3.680 

27 

21.250 

-6677.175 

-208.571 

-86.986 

0.233 

-6.294 

-3.679 

27 

23.750 

-6676.566 

-224.307 

-96.183 

0.253 

-6.294 

-3.679 

27 

26.250 

-6675.912 

-240.042 

-105.382 

0.273 

-6.293 

-3.679 

27 

28.750 

-6675.208 

-255.775 

-114.580 

0.293 

-6.292 

-3.678 
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XI 

Zteol  X  ' 

z 

X 

27 

31.250 

-SC74.447 

-271.501 

-123.775 

0.313 

-6.291 

-3.678 

27 

33.750 

•6C73.C41 

-287.226 

-132.969 

0.333 

-6.290 

-3.678 

27 

3C.250 

-6C73.785 

-302.953 

-142.162 

0.353 

-6.290 

-3.677 

27 

38.750 

-6471.874 

-318.676 

-151.352 

0.373 

-6.289 

-3.677 

27 

42.250 

-6670.519 

-340.684 

-164.219 

0.401 

-6.287 

-3.676 

27 

44.750 

-6669.493 

-356.401 

-173.406 

0.421 

-6.286 

-3.675 

27 

47.250 

-6668.419 

-372.117 

-182.595 

0.441 

-6.285 

-3.675 

27 

49.750 

-6667.290 

-387.827 

-191.780 

0.461 

-6.284 

-3.674 

27 

52.250 

-6666.105 

-403.533 

-200.963 

0.481 

-6.283 

-3.673 

27 

54.750 

-6664.890 

-419.244 

-210.150 

0.501 

-6.282 

-3.673 

27 

57.250 

-6663.604 

-434.943 

-219.329 

0.521 

-6.280 

-3.672 

27 

59.750 

-6662.278 

-450.645 

-228.506 

0.541 

-6.279 

-3.671 

28 

02.250 

-6660.909 

-466.343 

-237.688 

0.561 

-6.278 

-3.670 

28 

08.250 

-6657.398 

-503.999 

-259.701 

0.608 

-6.275 

-3.668 

28 

12.250 

-6654.903 

-529.091 

-274.373 

0.640 

-6.272 

-3.667 

28 

15.250 

-6652.943 

'547.904 

•28S.371 

0.6C4 

-6.270 

•3.666 

28 

17.750 

-6651.254 

•563.576 

-294.533 

0.484 

-6.269 

•3.665 

28 

20.250 

-6649.531 

-579.252 

-303.699 

0.704 

-6.267 

-3.664 

28 

22.750 

-6647.752 

•594.920 

-312.860 

0.724 

-6.265 

-3.663 

28 

25.250 

-6645.921 

-610.583 

-322.017 

0.744 

-6.264 

-3.662 

28 

27.750 

-6644.041 

-626.239 

-331.173 

0.764 

-6.262 

-3.661 

28 

30.250 

'6642.106 

-641.891 

-340.322 

0.783 

-6.260 

-3.660 

28 

32.250 

'6640.516 

-654.407 

-347.636 

0.799 

-6.259 

-3.659 

28 

35.250 

'6638.075 

-673.175 

-358.607 

0.823 

-6.256 

-3.657 

28 

37.750 

'6635.992 

-688.814 

-367.748 

0.843 

-6.254 

-3.656 

28 

40.250 

'6633.862 

-704.448 

-376.888 

0.863 

-6.252 

-3.655 

28 

42.250 

'6632.117 

-716.949 

-384.196 

0.879 

-6.251 

-3.654 

28 

44.750 

'6629.896 

-732.575 

-393.329 

0.899 

-6.249 

-3.652 

28 

48.250 

'6626.698 

-754.440 

-406.107 

0.926 

-6.246 

-3.651 

28 

50.750 

'6624.355 

-770.051 

-415.231 

0.946 

-6.244 

-3.649 

28 

53.250 

'6621.958 

-785.653 

-424.350 

0.966 

-6.241 

-3.648 

28 

55.750 

-6619.530 

-801.257 

-433.471 

0.985 

-6.239 

-3.647 

28 

58.250 

-6617.038 

-816.852 

-442.583 

1.005 

-6.237 

-3.645 
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X 

z 

39 

00.750 

-6614.496 

-832.438 

-451.693 

1.025 

-6.235 

-3.644 

39 

03.750 

-6613.430 

-844.904 

-458.979 

1.040 

-6.233 

-3.613 

39 

OS. 350 

-6606.585 

-879.167 

-479.003 

1.084 

-6.237 

-3.639 

39 

10.750 

-6603.849 

-894.732 

-488.099 

1.104 

-6.225 

-3.637 

39 

13.350 

-6601.073 

-910.294 

-497.193 

1.123 

-6.222 

-3.636 

39 

15.750 

-6598.338 

-925.844 

-506.282 

1.143 

-6.319 

-3.634 

39 

18.350 

-6595.336 

-941.380 

-515.359 

1.163 

-6.217 

-3.632 

39 

31.350 

-6591.818 

-960.028 

-526.254 

1.186 

-6.214 

-3.630 

39 

33.750 

-6588.830 

-975.561 

-535.328 

1.206 

-6.211 

-3.629 

39 

36.350 

-6585.793 

-991.084 

-544.398 

1.226 

-6.208 

-3.627 

39 

3S.750 

-6583.693 

-1006.595 

-553.459 

1.245 

-6.205 

-3.625 

39 

31.350 

-6579.563 

-1022.109 

-562.520 

1.265 

-6.202 

-3.623 

39 

33.350 

-6577.019 

-1034.511 

-569.768 

1.281 

-6.200 

-3.623 

39 

35.750 

-6573.783 

-1050.003 

-578.816 

1.300 

-6.197 

-3.620 

39 

38.350 

-6570.513 

-1065.497 

•587.865 

1.320 

-6.194 

-3.618 

39 

43.350 

•6565.169 

•1090.261 

-603.338 

1.351 

•6.189 

-3.615 

39 

44.350 

-6563.460 

-1102.640 

•609.560 

1.367 

-6.187 

-3.613 

39 

46.750 

-6559.037 

-1118.105 

-618.592 

1.386 

-6.183 

-3.611 

39 

49.350 

-6555.538 

-1133.561 

•637.619 

1.406 

-6.180 

•3.609 

39 

51.750 

-6551.994 

-1149.003 

-636.639 

1.425 

-6.177 

-3.607 

39 

54.750 

-6547.685 

-1167.531 

-647.456 

1.449 

-6.173 

-3.604 

39 

56.750 

-6544.768 

-1179.872 

-654.661 

1.464 

-6.170 

-3.603 

39 

59.350 

-6541.083 

-1195.293 

-663.665 

1.484 

-6.167 

-3.60C 

30 

01.350 

-6538.093 

-1207.625 

-670.860 

1.499 

-6.164 

-3.599 

30 

03.350 

-6535.081 

-1319.949 

-678.059 

1.515 

-6.161 

-3.597 

30 

05.750 

-6531. S66 

-1235.345 

-687.049 

1.534 

-6.158 

-3.595 

30 

11.350 

-6533.726 

-1269.198 

-706.810 

1.577 

-6.150 

-3.590 

30 

15.350 

-6516.351 

-1293.782 

-721.160 

1.608 

-6.144 

-3.586 

30 

17.750 

-6512.313 

-1309.144 

-730.124 

1.628 

-6.140 

-3.583 

30 

30.350 

-6508.208 

-1324.482 

-739.077 

1.647 

-6.137 

-3.581 

30 

23.350 

-6504.904 

-1336.754 

-746.237 

1.662 

-6.134 

-3.579 

30 

34.750 

-6500.726 

-1352.086 

-755.180 

1.682 

-6.130 

-3.576 

30 

37.750 

-6495.650 

-1370.469 

-765.907 

1.705 

-6.125 

-3.573 
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